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NOTATION 


A  a  vector  in  three  dimensional  space. 

A  an  operator  which,  when  operating  on  an  element  in  its 

domain,  gives  a  vector. 

ANx^=  a  column  matrix  with  dimension  Nxl  and  A^  as  its  qth 

element  (dimension  is  sometimes  omitted  when  understood) . 

A,.  »,=  [A  a  matrix  operator  with  dimension  M*N  and  A  as  its 

MxN  pq  MxN  pq 

element  in  the  pth  row  and  the  qth  column  (dimension  is 
sometimes  omitted  when  understood) . 


Superscript  *  denotes 
Superscript  T  denotes 


the  complex  conjugate  of  the  quantity, 
the  transpose  of  the  matrix. 
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I.  INTRODUCTION 

The  problem  of  penetration  of  electromagnetic  waves  through  an 
annular  aperture  in  an  infinite  conducting  screen  of  zero  thickness  is 
examined  in  the  frequency  domain.  The  formulation  includes  the  circular 
aperture  as  a  special  case.  The  generalized  network  formulation  for 
aperture  problems  by  Harrington  and  Mautz  [1]  is  used.  The  method  of 
solution  is,  in  general,  a  specialization  of  that  for  bodies  of  revolu¬ 
tion  by  the  above  authors  [2] . 

In  the  formulation,  the  equivalence  principle  [3]  is  used  to 
establish  an  integral  equation  for  the  unknown  magnetic  current  in  the 
aperture.  The  method  of  moments  [4]  is  used  to  solve  this  integral 
equation  with  expansion  functions  chosen  to  be  harmonic  in  <{>  (azimuth 
angle)  and  subsectional  in  p  (radial  variable  in  the  aperture) .  Because 
of  the  circular  symmetry,  Fourier  modes  are  decoupled  to  one  another 
and  only  matrices  of  considerably  smaller  sizes  need  :o  be  dealt  with 
one  at  a  time. 

Numerical  results  for  the  magnetic  currents,  radiation  patterns, 
and  transmission  coefficients  are  given  for  some  sample  cases  for  both 
normal  incidence  and  oblique  incidence. 

The  assumption  of  a  center  conductor  in  the  aperture  in  the 
present  problem  is  difficult  to  realize  in  practice.  However,  further 
investigations  of  problems,  such  as  wires  and  coaxial  lines  opening  into 
a  half  space,  may  utilize  some  of  the.  computer  programs  written.  Also, 
results  for  complementary  problems,  such  as  scattering  by  conducting 
washers,  can  be  obtained  from  the  solution. 
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II.  PROBLEM  FORMULATION 

The  problem  configuration  is  given  in  Fig.  1,  which  shows  an 
annular  aperture  in  a  perfectly  conducting  screen  with  and  RQut;  as 
the  inner  and  outer  radius  of  the  aperture,  respectively.  The  conduct¬ 
ing  screen  is  infinite  in  the  x  and  y  directions  and  has  zero  thickness. 
Impressed  sources  J1  and  M1,  which  produce  E!10  and  H10  in  the  absence 
of  the  screen,  exist  in  the  region  to  the  left  of  the  screen  (region  a, 

Z  >  0) .  The  region  to  the  right  of  the  screen  (region  b,  Z  <  0)  is 
source  free.  The  unit  normal  n  is  defined  on  the  x-y  plane  as  pointing 
away  from  region  a. 

The  equivalence  principle  is  utilized  to  separate  the  two  regions 
as  described  in  [1].  The  aperture  is  covered  by  a  perfect  electric  con¬ 
ductor  and  the  equivalent  magnetic  currents  are  used  to  produce  the  re¬ 
quired  tangential  electric  field  where  the  aperture  originally  existed. 

The  equivalent  situations  for  regions  a  and  b  are  shown  in  Fig.  2(a) 
and  (b) . 

For  region  a,  an  equivalent  magnetic  current  with  surface 
density  M  =  n  *  _E  (ji  denotes  total  electric  field),  which  is  nonzero 
only  over  the  aperture,  is  placed  just  to  the  left  of  the  plane.  This 
magnetic  current,  together  with  and  M*,  radiates  in  the  presence  of 
the  complete  conducting  screen  covering  the  entire  x-y  plane  to  produce 
the  total  field  in  region  a. 

For  region  b,  an  equivalent  magnetic  current  with  surface 
current  density-M  is  placed  just  to  the  right  of  the  plane  to  radiate, 
in  the  presence  of  the  complete  conducting  screen,  the  total  field  in 
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Fig.  2(a).  Equivalence  for  region  a. 
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Fig.  2(b).  Equivalence  for  region  b. 
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region  b.  The  use  of  H  in  region  a  and  -M  in  region  b  assures  the 
continuity  of  the  tangential  component  of  electric  field  across  the 
aperture. 

The  magnetic  field  in  region  a,  denoted  H  ,  is  the  sum  of 

that  due  to  the  impressed  source,  denoted  H*,  plus  that  due  to  the 

i  i 

equivalent  source  M,  denoted  H  (M) .  Note  that  H  is  produced  by  J 

i  a 

and  M  in  the  presence  of  the  conducting  plane  and  fl  ,  operating  on 

M,  gives  the  magnetic  field  in  region  a  due  to  M  in  the  presence  of  the 

conducting  plane.  Hence,  we  now  have  in  region  a, 

Ha  =  H1  +  Ha(M)  (1) 

The  magnetic  field  in  region  b,  denoted  Hb,  is  that  due  to 
the  magnetic  current -M,  denoted  Hb(-M).  Hence,  in  region  b,  we  have 

Hb  =  -  Hb(M)  (2) 

where  the  minus  sign  was  factored  out  due  to  the  linearity  of  the 
operator.  Again  note  that  -Hb(M)  gives  the  magnetic  field  produced 
by  -M  in  the  presence  of  the  conducting  plane  in  region  b. 

The  boundary  condition  remaining  to  be  satisfied  is  the  con¬ 
tinuity  of  the  tangential  magnetic  field  across  the  aperture.  This 
gives,  from  (1)  and  (2), 

-n  x  n  x  [fla(M)+flb(M)]  =  n  x  n  x  H1  (3) 

in  the  aperture  region. 

From  image  theory,  we  know  that 


I 
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n  x  H1  =  fi  x  2  H10 


(A) 


in  the  aperture  and  also  that  fla  =  -  2  fl  where  fl  is  defined  by 


H(M)  =  ( 


jwe 


— j  k | r-r ' 
r  )e  J  ' - 


M(r') 


J 

aperture 


r-r 


da' 


jwiJ 


[V^  •  M(j"')]  e 


-jk|r-r’ 


da'}  (5) 


aperture 


r  -  r' 


Here  k  is  the  wave  number,  £  and  y  are  free  space  permittivity  and 

o  o 

permeability,  respectively,  and  the  surface  divergence  operator  is 

V;  *  M(r’)  =  gfr  Mx(x,,y,)  +  g|rMy(x',y')  (6) 

where  M  and  M  denote  the  x  and  y  components  of  M. 
x  y  - 

From  (4)  and  (5),  (3)  is  now  written  as 


n  x  fi  x  H(m)  =  n  x  fi  x  |  h10 


(7) 


in  the  aperture.  Equation  (7),  together  with  the  definition  of  j} 
in  (5),  forms  the  basic  equation  for  determining  the  equivalent 
magnetic  current  M. 

The  method  of  moments  is  used  to  obtain  an  approximate 

solution  M  from  (7).  Let  M  be  a  linear  combination  of  a  finite 
— o  — o 

set  of  vector  "expansion"  functions,  defined  over  the  aper¬ 

ture  region 


M 
— o 


No 

y  v  m 

q-1  q_q 


(8) 


A  symmetric  product  of  two  vector  functions  over  the  aperture  is 


defined  as 
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<a,b>  = 


A  •  B  da 


(9) 


aperture 


A  finite  set  of  vector  "testing"  functions  {W  }  is  also  defined  over 

“P 

the  aperture  region.  We  assume  that  the  number  of  elements  of  {M^}  is 

N  ,  the  same  as  the  number  of  elements  of  {W  }. 
o  ~p 


A  matrix  equation 


ig»x»  v  - 1 

o  o 


(10) 


is  then  obtained  by  requiring  that  the  symmetric  product  of  each  with 

each  side  of  (7),  with  substituted  for  M,  be  the  same.  In  (10),  V 

is  a  column  matrix  with  dimension  N  whose  elements  are  coefficients  of 

o 

the  expansion  of  defined  in  (8),  that  is, 


V  '  'Vn*! 

o 


(11) 


The  I  is  also  a  column  matrix,  called  excitation  matrix,  of  dimen¬ 
sion  N 


where  I  is  defined  by 
P 


i  =  [i 


pJN  x  1 
r  o 


I  =  -<W  ,  hxfix-i-  H10> 
p  -P  -  -  2  - 


(12) 


(13) 


The  [Y  1  ..  is  a  square  matrix,  called  admittance  matrix,  of  dimen- 

pq  N  x  N 
o  o 

sion  N  *N  whose  elements  are  defined  by 
o  o 


Y  =<W,nxfixH(M)> 

pq  “P  ~  _  “<1 


(14) 


Furthermore,  if  the  testing  functions  are  all  tangential  to  the  x-y 
plane,  i.e. 

W  •  n  =  0  for  p=l,2,...,N 
— p  —  o 


(15) 
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then  (13)  and  (14)  can  be  reduced  to 


I  =  \  <W  ,  Hio> 
P  2  — p  - 


Y  =  -  <W  ,  H(M  )> 

pq  -p  -  ~q 

The  problem  now  is  to  determine  the  column  matrix  V  which, 
from  (10),  is  given  by 

V  =  [Y  ) _1X 

pq 

when  the  inverse  of  [Y  ],  denoted  [Y  ]  1,  is  assumed  to  exist. 

pq  pq 


III.  PROBLEM  SPECIALIZATION 


In  this  Section  the  expansion  and  testing  functions  used  in 
solving  the  problem,  and  also  the  type  of  incident  fields  under  con¬ 
sideration,  are  specified. 

(a)  Expansion  and  Testing  Functions 

The  expansion  and  testing  functions  are  described  in  polar  coordi¬ 
nates  in  the  x-y  plane.  The  domain  [Rj  ,  R  1  is  divided  into  M 

in  out 

equal  subintervals  with  a  uniform  length 


where 


A  -  d 

A  M 


d  =  R  -  R. 
out  in 


A  set  of  nodes  are  defined  at  the  boundaries  of  those  subintervals 


H  '  "in  *  “ 


for  Z  =  0, 1 , 2 , . . . ,M. 
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where 


For  expansion  functions,  we  use 


— q (P ’  «  =4  fsil(P)e 


q(s,Jl,  n)  =  (n+N)  (2M-1)  +  V^M-1)  +  l 
s  =  P,(f) 


i  =  1,2,. ...Lg 
n  =  -N,  -(N-l) , . . . , (N-l) ,N 


and  V  is  defined  as 
s 


if  s  =  p 


if  s  = 


L  is  defined  as 
s 


if  s  =  p 


if  s  =  <(> 


and  fg£(P)  is  defined  as 


V  i <0) 

fpK<P)  "  (p/d) 


V  Vp) 


where 


Vp)  ■ 


for  p.  ~  A  £  p  <  p. 


elsewhere 


p  <p)  -  p8(p-  f) 

2 

For  testing  functions,  we  use 

W  (p,<(>)  =  u^Ti(p)e‘ 
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where 

p(T,i,m)  =  (N+m)  (2M-1)  +  V^M-l)  +  1  (31) 

for  T  =  p,<l> 

i  =  1,2,... , 

m  =  -  N,  -(N-1),...,(N-1),N 


atid  8Ti(P)  is  defined  as 


gpi(p)  =  Clfpi(p) 


8*i(p)  =  C2  6(P"  (Pi  '  2}) 


where 


*1 


2 

Aird 


(32) 

(33) 

(34) 

(35) 


C2  Att 

and  6  is  the  Dirac  delta  function.  The  constants  c^  and  c^  are  used 
to  avoid  possible  unnecessary  common  factors  of  Y  and  T  in  their  final 
forms. 


(b)  Incident  Fields  ET*’0  and  H'*'0 

The  incident  fields  considered  here  are  plane  waves,  which  are 
good  approximations  for  fields  radiated  by  sources  distant  from  the 
aperture  region.  The  propagation  vector  k  is  assumed  to  lie  in  the 
x-z  plane,  and  the  angle  of  incidence,  0^,  is  defined  as: 

01  =  cos  1  (ufc  •  Gz)  0  _<  0*  £  j  (36) 

where 

Qk  =  k/Iisl  (37) 

For  each  propagation  vector  k,  two  unit  vectors  u  and  u  are  defined: 

-  -x  -// 
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u  =  u  x  u,  (39) 

-//  -y  -k 

The  Incident  plane  waves  considered  are  of  the  form 

^io  ^io  /\  -ik  •  r  rtrw 

E  =  E  u  e  J-  -  (40) 

—  — e 

/\ 

Hi0  =  x  Ei0  (41) 

no  " 

where  u  can  be  either  u  .  or  u  and  n  denotes  the  free  space 
— e  —  1  —  //  o 

impedance  /—  . 

Veo 

When  u  =  u  ,  the  incident  waves  above  are  said  to  be  of  vertical 
— e  —  i- 

(denoted  Ex)  polarization  and  the  polar  components  of  the  magnetic  field 
in  the  aperture  are: 

Hi°  =  -1^  cos  0i  cQS  ^  e-jkp  sin  61cos  <0  (42) 

P  no 

cos  0i  sin  ^  e-Jkp  sin  01  cos  4>  (43) 

*  no 

When  u  =  u  ,  the  incident  waves  in  (40)  and  (41)  are  said  to  be  of 
— e  — // 

horizontal  (denoted  E/y  )  polarization  and, 

Hi0  =  E^!  sin  0  e-jkp  sin  01  cos  0 
p  % 

Hf°  =  E^  cos  <{»  e"jkp  Sln  01  COS  ♦  (45) 

4»  nQ 

in  the  aperture  . 

Notice  that,  although  the  two  linear  polarizations  above  are  con¬ 
sidered  separately,  the  transmission  of  an  elliptically  polarized  wave 
can  be  obtained  from  the  results  of  these  separate  considerations. 


Appendix  A  gives  a  brief  discussion  of  this. 
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IV .  EVALUATION  OF  ADMITTANCE  MATRIX  AND  EXCITATION  MATRIX 


(a)  Admittance  Matrix 

In  the  preceding  sections,  the  general  formulation  of  the  problem 
and  the  specification  of  expansion  and  testing  functions  have  been  given. 
In  this  section,  approximations  used  in  the  computation  of  the  admittance 
matrix  [Y  ]  are  discussed.  These  approximations,  together  with  some  of 
the  intermediate  steps,  are  summarized  and  the  final  forms  of  the  matrix 
elements  are  given. 

(i)  The  surface  divergence  (6)  can  be  written  in  its  polar 


form  as 


A  =  —  -iL 

-  P  3P 


(PA  )  +  - 

p  p  H 


(46) 


where  A^,  A^  are  polar  components  of  A.  From  (22)  and  (46),  we  have 


*  *^(P'«  <t>')  =  Jn  P£(P')e^n<^  /p ’  for  s=4  (47) 


Recall  here  that  q  is  a  function  of  s,i  and  n  and  conversely,  s,£ 
and  n  are  each  functions  of  q,  as  shown  in  (23). 

To  compute  V't  •  M^<p '  , * )  for  s=p,  a  finite  difference  approxi¬ 
mation  is  used  to  spread  out  the  impulse  functions  encountered: 


V’ 

t 


*  2 


[P4(P’)  -  P£+1(p’)]ejn4,'/(Ap'/d) 


for  s=p 


(48) 


(ii)  The  admittance  matrix  elements  are  calculated  according 
to  (17)  because  the  testing  functions  (30)  are  tangential  to  the  x-y 
plane.  Equation  (17)  is  written  here  again  with  the  subscripts  of  the 
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for  n 


for 


and 


for 


and 


for 


and 


for 


Yn  = 


<A«J  IB"£] 

[°ul. 


=  -  N,  -(N-1),...,N  and  where 


<W  ,  .  .  ,  H(M  .  .  .  )> 

-p(P,i,n) - q(p,ft,n)  ' 


i  =  1,2,...,  M-l 
l  =  1,2,...,  M-l 


Biil  ~  '  %(p,i,n)’  -(^(<{),Jl,n))> 


i  =  1,2, ... ,  M-l 
SL  =  1,2 . M 


(55) 


(56) 


(57) 


<W  . .  .  , ,  H(M  .  „  J> 

—p(cf),i,n)  —  — q(p,  £,n) 


(58) 


i  =  1,2,.. . ,  M 
H  =  1,2,...,  M-l 


D.„  =  -  <W  , .  .  ,,  H(M  „  ,)> 

i£  -p((J),i,n)  — ~q  (4> ,  ^ ,  n) 


(59) 


1,2,...,  M 
1,2,...,  M 


The  column  matrices  Vn,s  and  In's  all  have  the  same  dimension 


(2M-1)  and  are  defined  as 


15 


16 


da 


out 


2n 


pdp 


dtp 


aperture 


aperture 


In 


0 


R  2n 

rp^f 


R. 

in 


0 


(64) 


(65) 


(v)  Inner  products  of  unit  vectors  are  summarized  here: 


A 

u 

A 

’  u  , 

=  u.  • 

u,  = 

COS  (<J)  *  -  4> ) 

HP 

-P' 

-4> 

“< t > 

A 

U  .  ' 

/V 

'  u  , 

=  -  u 

•  u.  , 

=  sin  (4>  *  — tp) 

-4> 

-p 

-p 

(66) 

(67) 


(vi)  It  can  be  shown  that,  for  all  admittance  matrix  ele¬ 
ments,  a  single  variable  ($'-<{))  in  the  angle  direction  is  needed 
during  the  integration  and  thus  a  single  variable  0,  instead  of 

(4> '  — <|))  ,  will  be  used  from  now  on. 

(vii)  New  variables  are  used  for  integrations  according  to 


p  =  At 

p'  =  At' 


and 


(9  +  1) 


(68) 

(69) 

(70) 


The  last  one  is  done  for  the  convenience  of  numerical  integration  by 
Gaussian  quadratures  in  the  range  of  <p  from  0  to  tt. 


(viii)  The  pulse  P  ^(p)  is  approximated  by  four  impulses 
during  the  integration. 


z+  y 


p„.  i<p)  4  v 

m=l 


■  ■ 


(71) 
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or,  equivalently 
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TT(f)  «  1  (79) 

where  A  is  the  wavelength  of  the  wave. 

The  final  forms  of  the  admittance  matrix  elements  are  now  given 
in  terms  of  the  submatrices  An,  Bn,  Cn,  Dn.  Note  that  although  approxi¬ 
mations  are  involved  as  described  in  (i),  (viii),  (ix)  and  (x),  equality 
signs  are  still  used  for  simplicity. 


.n  M 


pq  n 


=  ^T  K  Wm  {f  Fl(0m)[A  81  (p’q’£’m)] 


o  m=l 


&=1  f 


-tt  F  (9  )[g(p,q,6,m)  -  g(p,q,5,m)  +  g(p,q+l,5,m) 

J  K  J  Ill 


g(p,q+l,6,m)] } 


Bpq  =  -  TT  \  Wm  {f  W1  ^  h(p’q’*’m)]  +  t 

o  m=l  36=1 


F3(0mH8(p.q»6»ni)  -  g(p»q.5,m)]} 


c"  =  ~  l  w  {k(x  -  y)  F  (0  )  g  (p,q,5,m) 
pq  no  m=l  m  P  2  2  m  1 


K  F3(0mHg(P.q»5,m)  -  g (p , q+1 » 5 ,m)  ]  } 


Dpq  =  TT  \  Wm[jK(Xp  "  ¥  Fl(Gm)  h(p’q’5’m)  + 
rn  o  m=l 


Jk  F3^6m)  8(P’9»5’m)] 


where 


K  =  kA 
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0  's  and  w  's  (m  =  1,2, . ,N  )  are  the  abscissas  and  weight  factors 

mm  g 

of  Gaussian-quadrature  integration  of  order  as  described  in  (ix),  and 

F.(0  )  =  \  {cos [ (n+1)  f  (0  +1)]  +  cos [ (n-1)  £  (0  +1)]}  (85) 

1  m  1  i  m  i.  m 

F„ (0  )  =  4  {cos[  (n+1)  4  (9  +D1  ~  cos[  (n-1)  ^  +1)1)  (86) 

2  m  2  2  m  2  m 


and 


F_ (0  )  =  cos[n  y  (0+1)] 
J  m  2  m 


g(p,q,£,m) 


-jKR 


p£m 


dt 


q-i 


p£m 


(87) 


(88) 


^1  (p  >  q ,  £,m)  = 
2 


Xq+l/2  e“jKVm 


dt 


q-1/2 


p£m 


(89) 


where 


h(p,q,£,m)  = 


x  -j<R  „ 
q  te  pim 


dt 


q-i 


p£m 


(90) 


'P£m  =  >2  +  xp£  +  2t  xp£  Sin  [I  9m] 


x  =  x  +  q 
q  O 


q  =  0,  1,  y, . .  .  ,M 


(91) 

(92) 


Pp£  . 

V  =  A  =  Xo  +  p  + 


x  =  R.  /A 
o  in 


(£  -  2.5) 
4 


p  =  1,2,... ,M-1 
i  =  1,2, 3, 4 


To  evaluate  g,  g^  and  h,  an  approximation  of  e 
~2 

described  in  (x)  is  used  with  results  as  given  below: 


-j  KR 


p£m 


(93) 

(94) 


as 
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(b)  Excitation  Matrix  [I„l 

- p-s- 

The  excitation  matrix  is  considered  for  the  two  polarizations 
described  in  Sec.  III-(a). 


(i)  E ^-polar ization 

Substituting  (30),  (42)  and  (43)  into  (62),  one  obtains 


In  =  E10  ini 


pq 


.n 


pq 


1“  =  E10  inX 

0q  4>q 


(104) 

(105) 


where 


.ni  ,-ni  -1  ,-n+l  „i 

j.  =i  =  —  i  cos  0 

pq  pq  no 


q+l/2 


[Jn+1(Bt)  -  Jn_1(Bt)]dt 


q-1/2 


(106) 


C  =  ~  ^  TTo  j_n  Xq-l/2  COS  0l[Jn+l(PXq-l/2)+  Vl^q-l^  1 


B  =  <  sin  0 


(107) 

(108) 


The  integral  formula  for  Bessel  functions, 

2tt 
2tt 


.n 

J  (x)  =  J- 


e-Jn<t>  -  jx  cost})  d<J) 


(109) 


is  used. 


0 


It  can  be  shown  that,  when  j  <<  1  (typically  A  ~  ,  the  following 

approximation  for  inX  is  of  little  error  (typically  1%), 

pq 


i"1  =  TT  j‘n+1cos  91  [J  +1(Bx)  -  J  .(Bx  )] 
pq  p  n+1  q  n-1  q 


(110) 
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Now,  a  column  matrix 
-ni 


-nl 

v 


v 


-ni 


ni 

'pl 


ni 


M-l. 


ni 

'pl 


ni 

'<pM 


is  introduced  as  the  solution  of  the  equation 


§n  vni=  ini 


where 


rnl 

1 


Tni 

l 

P 


Tni 

ij. 


Lnf 

Pi 


.ni 

1pM-l 


.ni 


M>i 


ni 

4>M 


The  magnetic  current  in  this  case  can  now  be  written  as 

M  =  E  (m  u  +  m ,  u.) 

-o  P  -p  <f>  -1> 

where 


n  =  Tn  <P>  +  1  TL(P>COS  n<P 


P  P 
N 


n=l  PC 


1  =  l  T"  (P)sin  n(j) 


n=l 


(111) 


(112) 

(113) 


(114) 


(115) 


(116) 
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is  introduced  as  the  solution  of  the  equation 


^,n  -n//  -n// 

V  v  =  1 


(127) 


where 


-n//  -rn//  .n  u 

i  =  l  l  , 

P  Pi 


(128) 


.n  // 
1pM-l 


The  magnetic  current  in  this  case  can  be  written  as 


w  _io,  //  ~  .  ~  \ 

M  =  E  (m  u  +  m ,  u.) 
— o  p  —  <p  ~<p 


(129) 


where 


=  l  Tds(p)  sin  n<)) 

p  n=i  Ps 


(130) 


mJ  ’  TJ(P)  +  J1TJC(P)  COS  n<t> 

T”8(P)  '  ^  X  Vl/2(P)  ”-1>2" 
*P>  -  ”,  %  VP) 

q=l 


(131) 


(132) 


(133) 


T"  (P)  =  2  l  vn"  P  (p) 
4>c  <t>q  q 


n=l,2, . , 


(134) 


The  symmetry  used  in  this  case  is 


— n//  -n// 

V  =  -  v 

p  p 


(135) 


— n//  -n// 

v  *  =  v. 

<p  <P 


(136) 


V.  GAIN  PATTERN  AND  TRANSMISSION  COEFFICIENT 


Formulas  for  the  two  transmission  characteristics  are  given 
here  in  terms  of  the  magnetic  current  coefficients. 


(a)  Power  Gain  Pattern  on  the  Transmitted  Side 

The  power  gain  pattern  G(0,4>)  on  the  transmitted  side  is 
defined  as 

2 

G(9.<0)  =  Ttr2  C  |F  |  2  +  | Fq |  2 ]  /P 
o  ^ 

where  F. ,  F are  components  of  the  vector  potential  F, 


-jkr 


—  2vr 


M 
— o 


(r')e"jkr’COS  *  da- 


aperture 


in  the  far  zone  and 


(137) 


(138) 


cos  5  =  cos  9  cos  0'  +  sin  0  sin  0'  cos(<}>'  -  <J>)  (139) 


while  (r,0,<}>),  (r'  ,0'  ,$')  denote  the  field  point  and  source  point, 
respectively. 

P  is  the  time  average  of  the  total  power  transmitted  through 
the  aperture  and  is  given  by, 


I 


Re(n  •  E  x  H)  da 


or, 


aperture 


Re (M  •  H)  da 


aperture 


(140) 


(141) 


From  (2),  (8)  and  (141),  we  obtain 
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P  =  Re  l  V*  V  <M*,  H(M  )> 
ij  3  3 


(142) 


From  Section  III-(a),  we  notice  that  when  the  subsections  are  small, 

*  ^  e- 

in  (142)  can  be  approximated  by  with  a  constant  factor  and  an 
approximate  expression  for  is  obtained, 

.2  N  M-l  .  M  . 

P  =  Re  l  [  [  MV“  l“  +  [/  l“]}  (143) 

C  2  n-N  q-1  Pq  Pq  q-1  *q  ** 

From  (104),  (105),  (112)  and  (143),  we  obtain,  for  the 


Ej^-polarization, 


M  vni 

Pt  4  Re  {|Eio|W  l  en 

n=0  -ni 

v . 


(144) 


where  £  =  1  for  n  =  0  and  £  =  2  for  n  >  1.  From  (122),  (123),  (127) 

n  n  — 

and  (143),  we  obtain,  for  the  E  -polarization 


Pt  =  \  Re.{|Eio|2TiA2  l  £n 

n=0 


„  ~n// 
M  vp 


(145) 


Another  method  for  calculating  P  is  to  integrate  the  far  field 


radiation: 


>  .JeL[  - 

2n  J 


r2sin  0  d0  [  | F  |  2  +  | FQ |  2]d4> 


(146) 


The  ^-integration  in  Eq.  (146)  is  trivial  because  of  the  orthogonality 
of  the  trigonometric  functions.  The  0-integration  is  done  by  quadratures. 


F.,  Fq  used  in  (137)  and  (146)  are  given  in  terms  of  the  magnetic 

(p  U 

current  coefficients  in  two  separate  polarizations: 
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(ii)  E  -  polarization 

For  incident  fields  defined  by  (40)  and  (41)  with  u  =  u, , 

— e  — // 

substitute  (130)  to  (134)  into  (138),  assume  the  condition  of  (79),  and 
obtain 


F.  -  El0f? 


F  =  Elof" 

6  r0 


(151) 

(152) 


where 


2  -jkr  N  M-l 

f,  =  -  -  - -  y  e  cos  n^j^  7  vn//[J  ,  (kx  sin  0)  +  J  (kx  sin  0)] 

(J)  2r  n^Q  n  T  J  ^  pq 1  n+lv  q  n-1  q 


M 


+  jn+1  VJq[Xq-l/2(Jn+l(KXq-l/2Sln  6)  "  Jn-1 ' 9)) 


+  (2Jn(KXq-l/2Sin  0)  "  Jn+2(KXq-l/2sin  0) 


-  J  „(<x  /9sin  0))]}  (153) 

n-Z  q-l/z 

2  -jkr  N  M-l 

fg - ^ -  l  sin  n^-j1^  l  v^[Jn+1(  x  sin  0)  -  ^(^  sin  0)] 

n=l  q=l  t  “ 


M 


"  j  V<tq[Vl/2(Jn+l(<Xq-l/2Sln  6)  +  Jn-1 ' (lCXq-l/2Sin  6)) 


~  24°  (Jn+2(KXq-l/2Sln  0)  "  Jn-2(KXq-l/2Sin  0))]}  (154) 


(b)  Transmission  Coefficient 

The  transmission  coefficient  is  defined  as  the  ratio, 


TC  = 


(155) 


in 


where  P  is  defined  in  the  last  section  and  P.  is  the  time  average  of 
t  in 


.  ' —  — 

jr&ik  V 

'  . 


the  total  power  incident  upon  the  aperture  by  the  incident  sources  when 
the  incidence  is  normal.  It  is  given  by. 
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or 


I E'*'0 1  ^  2  2 

3.  =  - L-  u(R  -  R.  ) 

in  2n  out  in 

o 

I E^° I ^  2 

5,  =  ,  ‘  ~  M(M  +  2x  ) 

in  2n  o 


(156) 


(157) 


VI.  NUMERICAL  RESULTS  AND  DISCUSSION 

Numerical  results  for  a  number  of  cases  are  given  in  this  section. 

Several  examples  are  given  first  with  normal  incident  waves  where  only 

the  n=l  mode  is  needed  and  only  one  polarization  needs  to  be  considered. 

E^-polar ization  is  chosen  for  these  examples.  The  results  for  the  other 

polarization  are  easily  obtained  from  this  polarization.  Results  for 

oblique  incidences  are  illustrated  for  two  typical  apertures  (one 

small  and  one  intermediate).  Both  polarizations  are  considered  and 

an  adequate  number  of  modes  are  used. 

In  the  figures  showing  results  for  the  magnetic  currents, points 

are  plotted  according  to  the  actual  representation  at  the  centers  of  each 

subsection  (p-component  is  marked  by  symbol  A,  <}>-  component  is  marked 

by  symbol  CD) .  Note  that  subsections  for  the  two  polar  components  are 

shifted  from  each  other  by  Points  are  then  connected  by  straight  line 

segments.  At  p  =  R  ,  the  value  zero  is  plotted  for  the  p-component 

corresponding  to  the  boundary  condition.  At  p  =  R^,  except  for  the  case 

when  R.  =0  and  n  =  +  1,  the  value  of  the  p-component  should  also  go  to 
in 

zero.  This  is  discussed  in  previous  work  [5] .  For  the  exception  just 


mentioned,  no  value  is  plotted  at  p  =  R.  . 

in 

The  power  gain  patterns  are  computed  in  both  the  E-plane  and 
the  H-plane  at  elevation  angles  one  degree  apart.  The  transmission 
coefficient  is  computed  for  both  Ex  and  E  -  polarization  at  angles 
of  incidence  15°  apart. 

Figure  3  shows  both  components  of  the  magnetic  currents  for 

normal  incidence  on  an  aperture  with  R  =  .051,  while  R.  varies 

out  in 

from  0  to  .  01A  to  .  03A.  Figures  A  to  6  show  both  components  of  the 
magnetic  currents  for  circular  apertures  of  sizes  Rout  =  -25A,  .5A  and 
1.5A  respectively. 

These  results  (Figs.  4-6)  were  compared  with  previous  work  [7]. 

Better  convergence  is  observed  for  the  first  two  cases,  while  in  the 

third  case  (R  =  1.5A),  our  result  is  quite  different  from  theirs  and 
out 

seems  to  be  more  accurate. 

In  Figs.  7  to  11,  results  for  oblique  incidences  on  a  small 

aperture  (R  =  .02A,  R.  =0)  are  shown.  For  this  small  an  aperture, 
out  in 

the  n=l  mode  dominates  for  E^-incidence,  while  both  n=0  and  n=l  modes 
are  important  for  E//  -incidence.  This  result  can  be  seen  from  examining 
the  Bessel  functions  and  the  matrix  elements.  Both  components  of  n=l 
mode  are  shown  for  E1  -polarization  and  for  a  number  of  angles  of  inci¬ 
dence  in  Figs.  7,  8.  Figures  10,  11  are  the  equivalences  of  Figs.  7, 

8  for  E^  -  polar ization  while  Fig.  9  shows  the  circulating  current 
(n=0  mode)  for  this  polarization. 

Gain  patterns  for  this  aperture  are  also  computed.  For  E^- 
polarization,  a  simple  dipole  pattern  is  observed.  A  change  in  the 
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angle  of  incidence  only  changes  the  magnitude  of  the  dipole  (the  n=l 
mode)  and  doesn't  change  the  pattern  shape.  These  simple  dipole 
patterns  are  not  shown,  since  they  are  well  known.  For  E^-polar ization, 
since  both  dipoles  exist  (the  n=0  and  n=l  modes),  the  patterns  are  more 
complicated  and  are  shown  for  a  number  of  oblique  incidences  in  both 
planes  in  Figs.  18  and  19. 

In  Figs.  12  to  17,  magnetic  currents  for  oblique  incidences  on 

an  aperture  of  medium  size  (R  _  =  .  25A,  R.  =  0)  are  shown.  At  45° 

out  in 

incidence,  the  magnitudes  of  the  first  few  dominant  modes  are  compared 
for  both  components  and  for  both  polarizations.  The  same  is  done  for 
90°  incidence,  except  in  this  case  E^  -polarization  doesn't  excite  the 
aperture . 

Gain  patterns  for  oblique  incidences  on  this  aperture  are  shown 
in  Figs.  2 0  to  23,  again,  for  both  E  and  H-planes  and  for  both  polari¬ 
zations. 

The  transmission  coefficients  versus  angle  of  incidence  are 
plotted  in  Figs.  24  and  25  for  the  two  apertures,  both  polarizations 
considered.  For  the  case  where  R  =  .  02A,  the  result  compare  very 
well  with  formulas  by  Bethe  [6] .  A  slower  convergence  than  that  for  the 
currents  was  observed,  and  the  pulse  representation  of  the  edge  behavior 
of  the  ({(-component  of  the  current  is  considered  the  reason  for  this. 

Power  gain  patterns  shown  in  the  figures  are  normalized  with 
respect  to  the  maximum  gain  in  each  particular  plane.  Table  1  retains 
the  information  needed  when  these  patterns  are  normalized. 
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5.  First  mode,  p  and  <J>  components  of  magnetic  current,  |  Tpc  | 

and  ItL |  for  a  circular  aperture  of  R  =  . 5A, normal 
1  <j>s'  r  out 

incidence. 


RRD I US/R 


IX  2.0) 


Fig.  6.  First  mode,  p  and  components  of  magnetic  current,  |Tp  | 

and  It!  I  for  a  circular  aperture  of  R  =  1.5A,  normal 
'  As1  r  out 

inc idence. 


CURRENT 


o  .  1 


.9  1.0 


RflDIUS/R0UT 

g.  7.  First  mode,  p  component  of  magnetic  current  for  EA- 

polarization,  -jTj  for  a  circular  aperture  of  R  =  .02 
Pc  out 

angles  of  incidence  (a)  0°,  (b)  30°,  (c)  45°,  (d)  60°  and 

(e)  75°. 

(X  .35) 


g.  8.  First  mode,  4)  component  of  magnetic  current  for  EA~ 

polarization,  ,  for  a  circular  aperture  of  R  =  .02 
angles  of  incidence  (a)  0°,  (b)  30°,  (c)  45°,  (d?UEo°  and 
(e)  75°. 
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RRDIUS/Rout 

Fig.  9.  Zeroth  mode,  $  component  of  magnetic  current  for  E  - 

polarization,  -T°,  for  a  circular  aperture  of  R  =  .02A 
angles  of  incidence  (a)  30°,  (b)  45°,  (c)  60°,  ancf  (d)  90' 


(X  .  12) 


Fig.  10.  First  mode,  p  component  of  magnetic  current  for  E/(,  - 

polarization,  jT^  ,  for  a  circular  aperture  of  R  =  .01 
angles  of  incidence  (a)  30°,  (b)  45°,  (c)  60°  ancfU^d)  90' 
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Fig.  11. 


First  mode,  (p  component  of  magnetic  current  for 
polarization,  jl4  ,  for  a  circular  aperture  of  R 
angles  of  incidence  (a)  30°,  (b)  45°,  (c)  60°  an 


(X  1.51 


Fig.  12.  Zeroth,  first  and  second  mode,  p  component  of  magnetic 

current  for  EA -polarization  (a)  |T°|,(b)  lTpcl>  a°d 

(c)  |T2  I,  for  a  circular  aperture  of  R  =  .25X, 

1  pc  1  out 

45°  incidence. 
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Fig.  17.  Zeroth,  first  and  second  mode,  <p  con 
current  for  -polarization,  (a)  |l 
and  (c)  lT0cl>  for  a  circular  aperti 

90°  incidence. 


Curve 


ma 


-(c) 

-(d) 

-(e) 

-(f) 

19-(a) 

-(b) 

-(c) 

20 


G(0  ) 

max 


1.874 

2.195 

2.430 

2.581 

2.662 

2.687 

1.480 

1.341 

1.225 

2.331 


Curve 

G(8  ) 

max 

21-(a) 

2.358 

-(b) 

2.345 

-(c) 

2.331 

-(d) 

2.323 

22-(a) 

2.312 

-(b) 

2.191 

-(c) 

2.111 

2  3- (a) 

2.309 

-(b) 

2.184 

-(c) 

2.102 

Table  1.  Maximum  gains  in  each  normalized  gain  patterns  and  ele 
vation  angles  at  which  they  occur. 
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VIII.  CONCLUSION 


A  numerical  solution  for  the  problem  of  electromagnetic  trans¬ 
mission  through  an  annular  aperture  in  an  infinite  conducting  screen 
of  zero  thickness  is  developed.  The  moment  method  is  used  to  solve 
the  integral  equation  for  the  equivalent  magnetic  current  in  the  aper¬ 
ture.  The  expansion  functions  are  subsectional  in  p  and  harmonic  in  <)>. 
The  number  of  subsections  used  is  determined  by  the  size  of  the  aperture. 
The  number  of  Fourier  modes  used  can  be  determined,  in  general,  from 
the  size  of  the  aperture,  the  angle  of  incidence  and  an  inspection  on 
the  Bessel  functions.  The  results  are  observed  to  converge  with  re¬ 
spect  to  the  number  of  expansion  functions  and  some  of  the  results  are 
compared  with  previous  literature. 

An  attempt  was  made  to  use  triangle  expansion  functions  combined 
with  the  Fourier  modes  to  solve  the  same  problem,  however,  poorer  re¬ 
sults,  especially  for  small  apertures,  were  obtained.  The  edge  condi¬ 
tion  of  the  <J)-component  of  the  magnetic  current,  and  the  fact  that 
triangles  are  less  independent  from  each  other  (they  are  actually 
2-subsectional)  than  the  pulses  are  suspected  to  be  the  cause  of  the 


failure. 
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Appendix  A 

FORMULAS  FOR  ELLIPTICALLY  POLARIZED  INCIDENCE 


The  general  form  of  an  elliptically  polarized  plane  wave  inci¬ 
dence  can  be  written  as: 


_io  „  /  '-i  ~  x  -jk  •  r 

i  =  E0(el  u//  +  e2  ui)e - 

(A-l) 

where 

the  normalization  is  done  wuch  that 

|e1|2+  |e2|2  =  1 

(A-2) 

The 

equivalent  magnetic  current  for  this  incidence  is 

simply 

the  linear 

combination: 

M 

=  E0[(eimj  +  e2mJ)Gp  +  (e^J  + 

(A-3) 

—0 

where 

m  A, 

P 

m"  ,  mf,  mf  are  defined  in  Sec.  IV-b. 

P  <j>  0 

The 

far  field  components  are  linear  quantities  too: 

Fe  "  Eo(eif0  +  e2fp 

(A-4) 

%  -  w; + v  j> 

(A-5) 

where  f„  ,  f„,  f",  ,  ff  are  defined  in  Sec.  IV-b  also. 

0099 

The  transmission  coefficient  is  a  power  characteristic,  but. 


since  the  cross  products  don't  contribute  (which  can  be  seen  from  the 
symmetry),  the  total  transmission  is  still  a  linear  combination  of  those 
contributed  by  the  two  polarizations: 
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Re  J  Z  en  (ei 


T  T 

-n<r  *  ru  -nil*  ) 

/  M  v 

P  P  ! 

Til//  ,  TlU.I 

i  +  e2  1  )> 


M(M  +  2x  ) /r\ 
o  o 


(A-6) 


The  normalized  constants  e^  and  e^  for  an  elliptically  polarized 
plane  wave  can  be  obtained  from  the  ratio  of  the  length  of  the  minor 
axis  to  that  of  the  major  axis,  r,  the  angle  between  the  major  axis  and 
the  vector  ,  ip,  where  0  _<  iJj  _<  tt,  and  (when  the  wave  is  not  linearly 
polarized)  the  helicity  of  the  wave. 

An  elliptically  polarized  wave  is  said  to  be  of  positive 
(negative)  helicity  if  the  electric  field  vector  rotates  counterclockwise 
(clockwise)  when  the  observer  is  facing  into  the  oncoming  wave. 

If  r  ^  0  and  r  ^  1,  e^  and  e^  can  be  obtained  as: 


el  rl 


(A-7) 


e2  =  r2  e 


(A-8) 


where 


1.2,  .  cos  2\\) 

ri  =  /sln  *  +  7T^ 


1  +  r 


(A-9) 


r2  =  /cos  v 


2 ,  cos  2 
V  “  .  2 


1  +  r 


(A-10) 


-1  r sin  2'1'(1  -  r  ) 


cos  [ 


2rir2  (1  +  r  ) 


-1  rsin  2^(1  -  r  ), 
cos  [ - 3 - x-\ 

2rlr2  (1  +  r  ) 


]  -  tt  for  positive  helicity 


(A-ll) 


for  negative  helcity 


.  •  J'jjSSrA,  ’ %  V  ‘ ' 

' 
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I 

I  For  the  two  special  cases: 


(i)  When  r  =  1,  the  plane  wave  is  circularly  polarized  and  we  have. 


1 

e,  =  - 

1_ 

(A-12) 

1 

1 

/2 

1 

1  1 

1  zl 

for  positive  helicity 

| 

1 

1  /2 

(A-13) 

i 

for  negative  helcity 

1 

/ 

(ii)  When  r  =  0,  the  plane  wave  is  linearly  polarized  and  we  have. 
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Append ix  B 

COMPUTER  PROGRAM  LISTING  AND  DESCRIPTION 

Two  programs  are  given  in  this  section.  One  solves  for  the 
equivalent  magnetic  current  and  the  other  computes  the  power  gain  and 
transmission  coefficient  from  this  current.  Section  B-1.4  gives  the 
description,  listing  and  sample  output  of  the  first  program.  Sections 
B-l.l  to  B-1.3  give  the  description  and  listing  of  the  subprograms 
used.  Section  B-2.2  describes  and  lists  the  second  program  with  its 
sample  output.  The  subprogram  used  is  given  in  B-2.1. 

B-l.l.  Subroutine  CAD 

Subroutine  CAD(NP ,N1 ,N2 ,RD,RI ,LL,W,Q)  computes  submatrices 

Y  ,  Y  , . . . ,Y  of  the  admittance  matrix.  The  constant  factor  — 

o 

is  left  out  for  all  submatrices.  Elements  are  stored  in  Y  by  columns 

,-Nl 

and  one  submatrix  after  another,  starting  with  the  lowest  mode  (Y  ). 

The  input  variables  are  defined  in  terms  of  the  variables  introduced  in 
previous  sections  as 

NP  =  M,  number  of  subdomains  used. 

RO  =  R  ,  outer  radius  (in  wavelengths)  of  the  aperture, 
out 

RI  =  R.  ,  inner  radius  (in  wavelengths)  of  the  aperture, 
in 

LL  =  N  ,  order  of  the  Gaussian-quadrature  integration  in  <J>. 
g 

W(m)  =  w  ,  weight  factor  of  the  Gaussian-quadrature  integration  in  <)>. 
m 

Q(m)  =  0  ,  abscissa  of  the  Gaussian-quadrature  integration  in  <J>. 
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The  output  Y  is  transferred  to  the  main  program  in  the  common  block 
BLK2  by 

COMMON  /BLK2/Y 

to  avoid  unnecessary  storage  space. 

Minimum  allocations  are  given  by 
COMPLEX  Y (NY) 

DIMENSION  X(6) ,  Q(LL) ,  W(LL) ,  SN(LL) , 

CS(LL),  Ul(NL),  U2 (NL) ,  U3(NL),  D(NT) , 

T(NT) ,  Rl(NT),  R2 (NT) ,  R3 (NT) ,  R4 (NT) 

where 

NY  =  (N2  -  N1  +  1) (2M-1) 2 
NL  =  (N2  -  N1  +  1)LI. 

NT  =  6  LL 

The  matrix  elements  An  ,  Bn  ,  Cn  and  D°  are  computed  according 

pq  pq  pd  pq 

to  (80)  to  (83).  The  double  DO  loop  38  covers  the  ranges  of  both  sub¬ 
scripts  p,q.  DO  loop  36  does  the  Gaussian-quadrature  integration  in 

Ng 

4> ,  or,  equivalently,  the  summation  in  (80)  to  (83).  DO  loop  18 

4  m=l 

carries  out  the  summation  £  in  (80)  and  (81).  DO  loops  14,  17,  23, 

1=1 

28,  32  and  35  cover  the  range  n  =  Nl,  N1+1,...,N2  for  all  modes  con¬ 


sidered. 
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LISTING  OF  CAD 


SUBROUTINE  CADINP.Nl  . N2 . RO, R I , L L . W , Q ) 

COMMON  /BLK2/Y 

COMPLEX  Y (2500)  .C  .C J. CR ,C EX P , C E . E . F . G . H . A A, BB 

O I  MENS  ION  X(6>,0(30),SN(30».CS{30),Ut<360).U2( 360  > . U3 (360  >  . 
6  •  T  (  ISO  »,R1(IB0I.R2(180).R3<18  0),R4(180>,W(30) 

NR=NP— 1 
NN=NR*NP 
na=nn*  NR 

NB  =  NA«-NR 

NN2=NN**2 

JS=N1 ♦ 1 

JF=  N2  1 

JS1=JS+ I 

JF  1  =  JF  ♦  1 

JF2= JF  +  2 

JD= JF l- JS 

NT-NN2  *  JD 

P I =3 • 14159 

DD=RO-R  I 

A  K= ( 2*  P 14  DO ) / NP 

CJ=  <0. .  1  . J 

C=C J*AK 

X0=(R I *NP )/DD 

DO  1  1=1. LL 

Q  (  I)  =  Pl*(  Q(  1)41  )/2 

SN ( I)=SIN(Q(I  1  ) 

CS(  I ) =  COS ( Q ( I  I) 

K  =  0 

DO  2  I-l.LL 
DO  2  J  = 1. JD 
K  =  K*1 
U1 (K)=0. 

U2(K)=0. 

DO  5  1=1. LL 
K= JD* ( I-I ) 

DO  5  J= JS . JF2 
K  =  K  ♦  1 

AC  =  COS(  (  J  — 2  I  *Q{  It  ) 

IF< J.GT .JF ) GO  TO  3 
U1 (K)=U1 ( K) *AC/2. 

U2 I K )  =  U2( Kl-AC/2. 

IF<  J.EO.JSIGO  TO  4 
IF! J.GT  . JF 1  )GO  TO  4 
K 1 =K— 1 
U3IK1  )  =  AC 

IFIJ.LE • J  S 1 I  GO  TO  S 
K2=K— 2 

U1 CK2)=Ul C  K2) *AC/2. 

U2 (K2) =  U2 I K2»  *AC/2. 

CONTINUE 

DO  6  1=1.4 

X<  II  =  X0M  1-2.  5)  /4  . 

DO  7  1=1. NT 


01 180) 
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7  Y(  I  )=0  . 

00  38  1=1 .NP 
DO  8  K= 1 , 4 

8  XfK)— X(K)  +  1. 

X(5)=X0tl-.5 
X(6l=X0t| 4.5 
M=0 

DO  9  K=l.LL 
DO  9  L=  1.  6 
M  =  M41 

T ( M)  =  XO-X (L l*CS (K  ) 

D(M)=(X(L)*SN(K>)**2 
R3<  M)  =  SQRT<  T(  M>  **  24D(  M )  ) 

9  R41M)=SQRT(  {  T  (  M  )  4  .5  )  *  *2  4D  (  M  )) 

K2=NR4  ( 

K4=NB4I 

IF! I.EQ.NP1G0  TO  10 
K  1=1 
K3=NA4 I 

10  DO  38  J=l ,NP 
M=0 

DO  1 1  K=1  » LL 

DO  11  L=1.6 

M=M4l 

T ( M)  =  T ( M ) 4 1 • 

R  1  (MI  =  R3(  M) 

R2  (M)  =  R4J  M  ) 

R3(  M)  =  SORT  1  T<  M>  **24D<  M)  ) 

11  R4<M)  =  SQRT1 (T1M)4.5)**24D(M  )> 

K=0 

DO  36  M=1 .LL 
JLL=J0*(M-1 ) 

IF!  I  .EO.NPIGO  TO  19 

DO  18  MP= 1.4 

K=K41 

TT=ABS( T1 K )— . 5 ) 

T 1=TT4.5 
T2=TT- . 5 

TA=T14SQRT(Tl*»2*0(K) ) 

TB=T24SQRT(T2**2*D(  Kl  ) 

IF(T2.LT.O.)GO  TO  12 
HT=ALOG(TA/TB) 

GO  TO  13 

12  TB=  TB—  2*T  2 
HT=ALOG(TA*TB/D(K  )) 

13  CR=C*R2(K) 

CE=CEXP(-CR> 

H=CE*(  (  14CR)*<R3(K)-RICK)4X(MP)4CS<M»*HT)-C*(  J-.54X0)  ) 

JLB=JLL 

JNB=K3 

AA=-AK  *  W( M)*H/4. 

DO  14  JN=JS.JF 
JLB=JLB41 

Y(  JNBI  =  Y(  JNB  >  4  A  A*  U2  (  JLB  ) 

14  JNE= JNB4NN2 
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I F  C J.E0.NP1G0  TO  18 
SS=A0S<  T  (  K  1  ) 

S 1=  SS4  •  5 
S2=SS-  .5 

SA=Sl  ♦  SORT!  SI  **2+’0(  K 1  ) 

SB=  S2  ♦  SQR  T  C  S2  **  24  Ol  K I  » 

IF<  S2  .LT.  0  .  >GO  TO  15 
GT=ALOGC  SA/SB ) 

GO  TO  16 

15  SB=SB-2*S2 
GT=ALOG(SA*SB/D(K  >1 

16  CR=C*R3(K) 

CE=CEXP(-CR  ) 

G=CE*(  <  1  .♦CR)  *GT-C) 

JLB-JLL 

JNB=K 1 

AA=C*NP*W(M  >*G/4. 

00  17  JN=  JS  »JF 
JLB=JLB*t 

Y ( JNB)=Y( JNB1+AA4U1 (JL8) 

17  JNB=JN04NN2 

18  CONTINUE 
GO  TO  20 

1 9  K  =  K44 

20  K  =  K4l 

TT=ABS( T( K)-.5> 

T 1 =TT ♦ • 5 
T2=TT- .5 

TA=T1 +SORTC T1 **24D(K) ) 

TB=T2*SQRT<T2A*2«-0C  Kl  ) 

IF  (  T2  >LTaO>  )G.O  TO  2  1 
HT= ALOG ( T  A/TB ) 

GO  TO  22 

21  TB=TB-2*T2 
HT=ALOG(TA*TB/D(K >1 

22  CR=C*R2(K> 

CE=CEXP|-CR) 

H=CE*1  (  14CR)*(R3(  K  )— R  l(K)+X(5)  *CS(  M  1  *HT  )-C*  (  J-  .5+  XO  >  > 

G=CE*  (  I  I  .  4CR1  4HT- Cl 

JLB=JLL 

JNB=K4 

AA-C*W( Ml *< 1-.5+X0J *H 
B8=W< M ) *G/C 
00  23  JN=  JS »  JF 
JLB=JLB*1 

Y (  JNB  )  =  Y(  JNBM-AA4U1  (  JLB  »♦(  JN-  1  )  * ♦  2 *BB*U 3 (  JL B  ) 

23  JNB= JNB4NN2 
TT=ABS(TC  K  )  ♦  •  5 1 
T 1=TT 4.5 

T2=  T  T-  .  5 

TA=Tl ♦SORT ( Tl *424 01 K> I 
TB=T24SQRTC  T2**2*0l  Kl  1 
IFIT2.LT. 0. 1G0  TO  24 
GT  =  ALGG(T  A/TB > 

GO  TO  25 

24  TB=TB-2*T2 


GT=ALOGCT  A*TB/Df X ) ) 

25  CR=C*R4CK> 

CE=CEXPC-CR ) 

E=CE*C  <  l.«-CR»  *GT-C) 

IFC J.EQ.NP1G0  TO  29 
SS=ABS<  TC  K  1) 

Sl=SS+.5 

S2=SS- . 5 

SA  =  S1  ♦SQR  T ( SI *#2*DCK)  » 

SB  =  S2+S0RTC  S2**2*OC  K)  ) 

IFCS2.LT. O.IGO  TO  26 
GT=ALOGC  SA/SBI 
GO  TO  27 

26  SB-=  SB—  S2*  2 

GT= ALOG ( S A* SB/O ( K  )) 

27  CR=C*R3CK) 

CE=CEXPC-CR) 

F-=CEA  {  (  I  .  *CRi  *G  T—  C ) 

JLB-JLL 

JNB=K2 

AA=NP*W  (M  )  *AK*(  l-.5*XO)*F 
B8  =— NP*  M ( M) *< G-E> /AK 
00  28  JN=JS.JF 
JLB=JLB«-l 

Y( JNB )  =  VC  JNB) +AA*  U2  C  JLB >  ♦  C  J  N-  1  )  *  BB  *  J  3  C  J  LB  ) 

28  J NB=JNB*NN2 

29  K=K+1 

I F (  I.EQ.NP1GO  TO  36 
TT  =  ABS C T<  K)-. 5) 

Tl=TT«-.5 
T2  =  TT-  .5 

TA=Tl«-SQRT<  T1  **2+0C  K)  ) 

T8=T2«-SQRTC  T2**24  0C  K>  ) 

IFCT2.LT.  O.IGO  TO  30 
HT=ALOGCT A/TB) 

GO  TO  31 

30  TB=TB-2*T2 

HT  =  ALOG  CT  AATB/DIK  )) 

31  CR=C*R2(K> 

CE=CEXPC-CRI 

F  =  CE*C  C  1 . ♦CR) *HT— C) 

JLB=  JLL 
JNB=K3 

AA=-WC  M )*  CF-G )/AK 
OO  32  JN=JS.JF 
JLB  =  JLB  ♦I 

VC  JNB  )  =YC  JNB)  *C  JN-l  )  *  AA4U3C  JLB1 

32  JNU=JNB*NN2 
1FCJ.EQ.NP) GO  TO  36 
TT=ABS ( TC K  ) ♦•SI 
Tl=TT*.S 

T2=TT- . 5 

TA=ri«-SQRT(Tl*4240CK|  I 
TB=T2*SQRT«  T2  *4  24  01  K)  I 
IFCT2.LT. O.IGO  TO  33 
GT=ALOG(T  A/TB  ) 


m, " 


33  TB=T8-2*T2 

GT  =  ALOG  (  T A4TB/D (K > ) 

3*  CR=C*R4(K) 

CE-=CEXP  I— CR  ) 

H=CE4<  I  1  • ♦CR) *GT— C) 

JLB= JLL 
JNB-K1 

AA=-NP»W(  M)*<F-G-H*E)/C 
DO  35  JN=JS.JF 
JLB=JLBM 

Y<  JNB»  =  Y(  JNB)  *AA*U3(  JLB  ) 
35  JNB=JNB«-NN2 


36 

CONTINUE 

K4=K4+  NN 

I F  C  1  .EQ.NP)  GO 

TO 

37 

37 

K3=K3+NN 
IF(J.GE.NR)  GO 

TO 

38 

K2=K2 * NN 

IF ( I.EQ.NPJGO 

TO 

38 

38 

K 1— Kt f NN 

CONTINUE 

RETURN 

ENO 


a®5 
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B-1.2.  Subroutine  PEX 

Subroutine  PEX(AI ,RO,RI ,M,N)  computes  the  column  matrices 
i0i,  i11,...,  iNX  for  the  E^polar  ization  and  I0//,  i1^, . . .  ,IN//  for 
the  E^-polarization.  The  constant  factor  is  left  out  for  all 
matrices.  The  ?ni's  are  stored,  starting  from  I  ,  in  W  and  i  's 

in  Q.  The  input  variables  are  defined  as 

AI  =  01  ,  the  angle  of  incidence,  in  degrees 

RO  =  RQut_,the  outer  radius  (in  wavelengths)  of  the  aperture. 

RI  =  R  ,the  inner  radius  (in  wavelengths)  of  the  aperture. 

M  =  M  ,  number  of  subdomains  used. 

The  outputs  W  and  Q  are  transferred  to  the  main  program  in  the  common 
block  BLK1  by 

Common  /BLK1/W,Q 

to  avoid  unnecessary  storage  space. 

Minimum  allocations  are  given  by 
COMPLEX  W(NI) ,  Q(NI) 

DIMENSION  X(M-l) ,  Y (M) 

COMMON  /BLK4/BJ (N+2) 

where 

NI  =  (2M-1) (N+l) 

The  matrix  elements  of  Ini's  and  in//,s  are  computed  according  to 
(106),  (107),  (124)  and  (125).  DO  loop  9  adds  Jn(x^),  n=0, 1 , . . . ,N+1 , 
q  =  1,2,...,M-1,  one  at  a  time,  to  the  proper  elements  of  W,Q,  with  the 
proper  constant  factors.  DO  loop  14  has  the  same  function  for  Jn(xq  l/2^'s* 
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If  AI  =  0  (0^  =  0),  only  i^and  i^are  stored  in  W  and  Q  and  they 
are  stored  as  W(2M)  to  W(4M-2)  and  Q(2M)  to  Q(4M-2).  Simple  computation 
are  done  without  computing  the  Bessel  functions.  This  is  done  by  the 
control  statement 

IF (AI .  EQ.  0.)  GO  TO  15 


<0  a 


LISTING  OF  PEX 


SUBROUTINE  pex< M .ro.ri ,M,N) 

COMMON  /BLKt/W.  0/BLK4/B  JI  100) 

COMPLEX  WI300  )  . QI 30 0 ) . C J. CK .C A . T W » TQ 
DIMENSION  X  I  2  4  )  »Y  (  25  )  *• 

MM=M—  1 

MO=M*MM 

NP-N*1 

NP2=NP+l 

MN=MD*NP 

D=RO— R  I 

CJ=<0. . 1. ) 

XO=M*R I/O 
E  =  X  0—  .  5 
DO  t  I  ■=  1  *  M 
E=E«-l . 

1  VI  I  )=E 
E=XO 

DO  2  1=1. MM 
E=E«-l . 

2  XII)=E 

IFIAI • EQ. 0* )  GO  TO  IS 
PI=3. 141593 
A 1  =  AI*P  1/ 1  BO. 

A=COS ( All 

B=2.*PI*D*SINIA1)/M 
CA=CJ* A 
DO  3  1=1. MN 
U  I  I  )=  0  . 

3  Q I  I  I =0  . 

DO  9  1=1. MM 
T=B*XI  I  ) 

CALL  BES(  NP.T  ) 

K=  I 

IFIAI • E Q. 90 • ) GO  TO  4 
WIK)=-2.*CA*BJI2) 

4  CK= I . 

DO  9  J  =  1  • NP2 

KsK+HO 

CK-CK/CJ 

IFIAI. EQ. 90. >60  TO  5 
TW=CA*6JI J ) /C K 

5  TQ=-BJIJ)/CK 

IF (J.GT.N) GO  TO  7 
IFIAI .EQ. 90.) GO  TO  6 
M I K )=W I K) +  T» 

6  QIK|=QIK)-TQ 

7  IFIJ.LT .3 ) GO  TO  9 
L*K— 2*  MO 

IFIAI .EQ. 90.) GO  TO  8 
M(L)  =  W<  L) 4TM 
QIL)=QIL) *ro 
CONTINUE 


DO  14  1=1 ,M 

T=B*Y< I > 

CAUL  BES(NP.T) 

CK=1  . 

K= I4MM 

Q<K)=-2.*CJ*Yf l >*BJ<2) 

DO  14  J= 1 • NP2 

K=K+MD 

CK=CK/CJ 

IF  <A1  .EQ.  90.)  GO  TO  10 
T  W=— A*  Y ( I  ) 4BJ ( J l/CK 

10  TQ=CJ*Y(I )*BJ< J)/CK 
(F(J.GT.N)GO  TO  12 
(FIAI .EQ.90.IGO  TC  II 
W<K)=M(K)*T* 

11  Q(K)=Q(K)+TQ 

12  IFIJ.LT.3) GO  TO  14 
L=K— 2  AMD 

I F ( A  I .EQ. 90.) GO  TO  13 
«(L)=M(L)-TM 

13  Q(L)=QU.)+TO 

14  CONTINUE 
GO  TO  18 

15  K=*D 

DO  16  I =1 .MM 

K=K+1 

Q( K )=— C J 

16  W( K  )=1  . 

OO  17  1=1  ,M 

K=K41 
0<K)=Y(  I) 

17  «<K»=CJ*Y(I) 

18  RETURN 
END 
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B-1.3.  Subroutines  DECOMP,  SOLVE  and  BES 


The  subroutines  DECOMP  and  SOLVE  uses  the  method  of  Gaussian 
elimination  and  LU  decomposition  to  solve  a  linear  system  of  equations 
with  complex  coefficients.  These  subroutines  are  described  in  [8].  Only 
input  instructions  and  output  descriptions  are  summarized  here. 

The  input  into  DECOMP  (N,  IPS,  UL)  is  N  and  the  matrix  A  of  coef¬ 
ficients  of  the  set 


A  x  =  b 

of  N  linear  equations  stored  by  columns  in  UL.  N  >.  2.  The  output  from 
DECOMP  is  IPS  and  UL.  DECOMP  does  not  change  N.  SOLVE  (N,  IPS,  UL,  B,  X) 
uses  N,  IPS,  UL,  and  the  elements  of  b  stored  in  B  to  calculate  and  stored 
in  X  the  elements  of  the  solution  x.  SOLVE  does  not  change  any  of  the 
input  variables. 

Minimum  allocations  are  given  by 
COMPLEX  UL(N2) 

DIMENSION  SCL(N) ,  IPS(N) 

in  DECOMP  and  by 

COMPLEX  UL(N2),  B(N) ,  X(N) 

DIMENSION  IPS (N) 


in  SOLVE. 

Subroutine  BES(N,X)  computes  and  store  in  BJ  the  values  of 
Jq(x),  J1(x),  ...,  JN(x)  by  the  numerical  method  in  9.12  on  page  385 
of  [9]. 


Minimum  allocation  is  given  by 


COMMON  /BI.K4/BJ (100) 


— - TT 


LISTING  OF  DECOMP 


subroutinE  de co MP ( N , 1  P  S  »  UL  ^ 

COMPLEX  ULI2500), PIVOT, EM 
DIMENSION  SCL (50) .IPSI50I 
DO  5  1 =1 , N 
IPSC I >=I 
RN*0  • 

Jl-I 

DO  2  J=1,N 

ULM=ABSIREALI UL ( J I) ) ) ♦ A BS ( A  IM AG < UL I J t  I )  I 
JI=JI  «-N 

IFIRN-ULM)  1,2,2 

1  RN=ULM 

2  CONTINUE 
SCL I  I )  =  l,/RN 

5  CONTINUE 

NM l=N— t 
K2=0 

DO  17  K  =  l ,NMl 
a  I  g=o . 

DO  II  l— K • N 
IP=IPSI  I  ) 

IPK=IP  + K2 

SIZE=I ABSIREALIULI  IPK)  )  M-ABSI  A I  MAG  ( UL  (  I  PK  I  >  )  »  *SCL  I  I  P I 
IF  I S I ZE— B I G  I  11,11*  10 

10  B  IG  =  S1 ZE 
IPV=1 

11  CONTINUE 

IF  I  IPV— Kl  14,15.14 

14  J— I PS (  K  I 

IPSIK  )=IPS(  IPV) 

IPS!  IPV  >=  J 

15  KPP=I PSIK )+K2 
P I VOT=  ULI KPP) 

KP1=K«-1 

DO  16  I=KP1,N 
KP— KPP 
IP-=  IPSI  I  )  4X2 
EM=— UL I  IP l/PI VOT 
18  ULI  IP  I--EM 

00  16  J-KP l ,N 
I P— I P* N 

kp*kpfn 

ULI  I  Pt  —  UL  I  I  PI  ♦EM4ULIKPI 

16  CONTINUE 
K2=K2 ♦ N 

17  CONTINUE 
RETURN 
ENO 


*  u 


LISTING  OF  SOLVE 


SUBROUTINE  SOLVE!  N.  IPS  .  UL  .  0  .  <  I 
COMPLEX  UL  (  25  00  )«0(50)«X(50)»  SUM 
DIMENSION  IPSISO) 

NP1=N4 1 
IP=IPS( 1> 

XI  I )  =  B(  IP  ) 

DO  2  1=2. N 
IP= IPS! I > 

I  PB=  I P 
IM1=I- I 
SUM=0. 

DO  1  Js|,  IM1 
SUM=SUM*-UL  I  IP  I*  XI  J) 

|  I P= I P+  N 

2  XI  I  i  =  B  (  IPB  )  — SUM 

K2=N4< N-l I 
IP=IPS<  N) 4K2 
XI N)  =  X INI/ULI  IP  I 
DO  4  IBACK=2.N 
1=  NP 1 —  I  BACK 
K2=K2—  N 
IPI=(PS<I  >*K2 
IP  1=1  ♦  l 
SUM=0. 

IP=IPI 

DO  3  J= IP  1  *  N 
IP=IP4N 

SUM=SUM«-UL(  I P  I  ♦  X  (  J) 

XI  I >  =  t XI  I  )—  SUM I /UL I  IP  I  ) 

RETURN 
END 


LISTING  OF  BES 


SUBROUTlNE  BESIN.X* 

COMMON  /BLK4/BJI  l  0O)/BLK3/T0(  33  > 
P 1=3. 141593 
P I  2  =  2 • /PI 
P  I  4=P  I  /4  • 

PI 7=. 7 S»P  I 
MZ=10*N*2*IFI X< XI 
BJIMZ4-1 1=0. 

BJ1MZ 1= l.E-60 
ALP=0. 

IF IMZ-I MZ/2142)  13.14*13 

14  JT=— 1 

GO  TO  15 
13  JT= 1 

15  M2=MZ— 2 


X2  — 2  • / X 

DO  16  K  =  1 . M2 

MK=M2-K 

B J<  MK  »  =  X2 *FLO AT  <MK) *BJ<  MK  M  1-BJ ( MK+2 I 

JT=-JT 

S=  1+JT 

ALP=ALP*B J<  MK)*S 

16  CONTINUE 

BJ(ll  =  X2*aJI2»-BJ(3) 

AL  P=ALP*B J (  l  » 

NP-N-M 

DO  17  K  =  1  ,NP 
BJ(K)=BJ( K ) /ALP 

17  CONTINUE 
RETURN 
END 
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B-1.4.  Main  Program  1. 

This  program  treats  an  annular  aperture  with  a  number  of  inci¬ 
dences.  For  each  incidence,  the  admittance  matrix  is  the  same,  only 
excitations  vary.  Punched  data  cards  are  read  early  in  the  main  pro¬ 
gram  according  to  input  statements  in  the  following  sequence. 


101 

FORMAT 

(13) 

102 

FORMAT 

(2E14.7) 

103 

FORMAT 

(16F5.1) 

104 

FORMAT 

(2613) 

114 

FORMAT 

(10F8.5) 

READ 

(1 , 101)M 

READ 

(1,101) ICASE 

READ 

(1,101)LL 

READ 

(1,102)  R0,RI 

READ 

(1,103)  (AI(J), 

J=l, 

ICASE) 

READ 

(1,104)  (NS ( J) , 

J=l, 

ICASE) 

READ 

(1,104)  (NF ( J) , 

J-l. 

ICASE) 

READ 

(1,114)  (WW(J), 

J-l. 

LL) 

READ 

(1,114)  (QQ(J), 

J-l, 

LL) 

The  above  input  data  is  defined  as 

M  =  number  of  subdomains  used. 

ICASE  =  number  of  incidences  considered. 

LL  =  order  of  the  Gaussian-quadrature  integration 
in  for  admittance  matrix  elements. 

RO,RI  =  outer  and  inner  radius  (in  wavelengths)  of  the 


aperture,  R 


64 


AI(i)  =  angle  of  incidence  of  the  ith  incidence 
(in  degrees) . 

NS(i)  =  order  of  the  lowest  mode  used  for  the 
ith  incidence. 

NF(i)  =  order  of  the  highest  mode  used  for  the 
ith  incidence. 

WW(i)  =  weight  factor  of  the  Gaussian-quadrature 
integration  in  4). 

QQ(i)  =  abscissa  of  the  Gaussian-quadrature 
integration  in  <p. 

Minimum  allocations  are  given  by 

COMPLEX  Y (MY) ,  Z(MZ),  W(NI) ,  Q(NI) , 

PI (2M-1) ,  P2(2M-1),  V1(2M-1),  V2(2M-1). 

DIMENSION  AI(ICASE),  NS(ICASE),  NF(ICASE) 

IPS (2M-1) ,  IPSL(NI),  WW(LL) 

QQ(LL) 

where  MY  =  (max  (NF(i) )  -  min  (NS(i))  +  1)  (2M-1)2 

MZ  =  (2M-1) 2 

NI  =  (max  (NF(i) )  +  1)  (2M-1) 

DO  loop  3  determines  the  order  of  the  highest  and  the  lowest  modes, 
MA  and  MI,  needed  for  all  incidences  considered.  These  two  integers  are 
then  fed  into  the  subroutine  CAD  as  the  input  parameters  Nl  and  N2.  In  DO 
loop  6,  admittance  matrices  Y^,  Y^"*"^ , . .  . ,  Y^  are  taken  out  of  Y,  one  at 
a  time,  stored  in  Z,  and  fed  into  DECOMP,  the  output  Z  is  stored  back  into 
Y  and  the  output  IPS  is  stored  into  IPSL  for  each  mode.  This  is  done  early 
because  the  LU  decomposition  needs  to  be  done  only  once  for  all  excitations. 
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Different  modes  of  different  incidences  are  then  treated  in  the  double 
DO  loop  14  where  subroutine  PEX  is  called  to  generate  excitations  for 
both  polarizations  according  to  the  modes  specified  for  each  incidence. 
Different  modes  are  then  solved  for  each  incidence.  Here,  submatrices  for 
each  mode  are  again  taken  out  of  Y  and  IPSL  and  fed  into  SOLVE  together 
with  the  proper  excitation  matrices.  The  output  VI  is  the  solution  vnJ" 
and  the  output  V2  is  the  solution  vn//. 

Notice  that  the  current  solutions  are  printed  out  at  each  stage 
but  not  kept  in  a  permanent  storage.  Data  should  be  punched  out,  if 
needed,  inside  the  DO  loops. 
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PROGRAM  1 

COMMON  /BLK1/W. Q/ BL K2/ Y /BLK  3/ Y 0  <  33) 

COMPLEX  Y (250  0)  ,Z  !  2 50 0  )  .W  (  3 00  )  .  0 ( 300)  * P  t ( 50  )  »P2 ( 50 )  . V 1 ( 50  )  * V2 ( 50 )  . 
6CJ. VT 

DIMENSION  A  1(20). NS (20 ) .NF<  20 )«  IPS (50) .  I PSL ( 300 ) • WW ( 3 0) .00(30) 

10)  FORMAT (13) 

102  FORMAT  (  2E  1  4. 7) 

1  03  FORMAT  (  I6F5.1  ) 

104  FORMAT (26 13) 

105  FORMAT  < • 1  •  »  IX  .  •  ANNULAR  APERTURE  PROBLEM*.///) 

106  FORMAT (5X .'OUTER  RADI  US /WAVELENGTH *. 4X .* I NNER  R AD ( U S/ MA VELENG TH* . 4 
6X. 'NUMBER  OF  SU BOOM A  1 NS • , / . 1 1 X . E 1 2 .5 . 1 6 X, E 1 2 . 5 .  1 7X .  1 2  ,// ) 

107  FORMAT ( 5X , 'CASE'. 3X . • ANGLE  OF  I NC I OENCE ' * 3X . • MOOES  USED • , /. 16X.  •  ( O 
6EGREES ) • .// ) 

108  FORMAT! 5X . I  3.  10X.F5. 1 . 10X. I  3. *  TO'. 12) 

109  FORMAT! *1 • . 5X . 'EQUIVALENT  MAGNETIC  CURRENT:'.///) 

110  FORMAT !6X.' CASE',  13.//) 

111  FORMAT (/. 7X .' N= • ,  12 .//, 1 l X, 'VERT ICAL  POLARIZATION'./) 

112  FORMAT ( 1 1 X. 12 .5E1 2.4) 

1 13  FORMAT ! 1 1 X •  12 .2 X.  '0.'  .lOX.'O.'.lOX.'O.'  ) 

1 14  FORMAT! 10F8.5) 

115  FORMAT (/. 1 1 X, 'HOR IZCNTAL  POLARIZATION'./) 

116  FORMAT  ( ////.5X,  13  .'  POINTS  GAUSSIAN  QUADRATURE  INTEGRATION  IS  USED 
6*  ) 

REAO! 1  .  10  1  )  M 
READ!  1*10  1)  ICASE 
READ!  1  .  10  1  ILL 
READ!  1.102)R0.R( 

REAO! 1 .103)(Al(J).J  =  l , (CASE) 

REAO!  1  .  10  4)  (NS<  J)  ,J*1  .ICASE  ) 

REAO! 1 • 10  4) (NF! J)  •  J= 1 . 1  CASE  ) 

REAO! 1 . 1 l 4) ( WW( J) .J-l ,LL) 

READ!  1  *  11 4) (00! J)  .J-l ,LL) 

WRITE! 3.1 05) 

WRITE! 3.1 06)RO.RI .M 
W  R  I  TE  (  3 , 1  0  7  ) 

MM=M—  1 
MD=M+MM 
MD2=MD*42 
C J=(0 . *  1  .  ) 

X0=M*R I /( RO-R  1  ) 

00  1  J=  1.  ICASE 

1  WRITE  (3.1  08)  J.AK  J)  .NS!  J)  «NF(  J) 

WRITE! 3.1  1 6  ILL 

MA  =  0 

M  I  =  100000 
00  3  J*l.  ICASE 
IF(NS(J).  GE  .MI)  GO  TO  2 
M  I— NS  (  J  ) 

2  IF! NF ( J).LE.MA) GO  TO  3 
MA*NF(  J) 

3  CONTINUE 

CALL  CAD(  M.MI  .MA.  RO.R  I.LL  •  WW.  QQ  ) 

K=  0 


’  v  ■'  ■ 


4 


KK=0 
L  =  0 

ND=MA—  M  I  1 
DO  6  J=  l » NO 
00  4  N=l,MD2 
K  =  K  + 1 
Z!N)=Y(K) 

CALL  DECOMPIMD.  IPS.  Z) 

DO  5  N* », MD2 
KK  =  KK*  1 

5  Y(KK)=Z(N) 

OO  6  N= 1 , MO 
L=L*I 

6  I  PSL(  L  )=l  PS  IN  ) 

WRITE! 3.1 09) 

00  14  J  =  l  . I  CASE 
WRITE! 3.1  10) J 
JS=NS( J) 

JSl=.JS*l 
JF=NF( J) 

JF  l  =  JFtl 

AN=  A I ( J ) 

CALL  PEX( AN.RO.RI  .M. JF) 
K=MD2*<  JS-Mt) 

L=MO* JS 
LL=MD*  C  JS-M I ) 

OO  14  NsJSl.JFl 
N  l=N-  1 

DO  7  1=1,  MD2 
K  =  K+l 

7  Z I  I )=V (  K) 

00  8  1=1. MO 
L=L41 
LL=LL» 1 

PI  (  I  )  =  W(L  ) 

P2< I )  =  0 (L ) 

8  IPSII  )  =  IP  SL ILL) 

CALL  S0LVE1MD.IPS.Z.PI . VI ) 
CALL  SOLVE  I  MO • ( PS ,Z.P2. V2 ) 
WRITEI3.1 11)N1 
OO  9  I =1, MM 
VT  =  VI  I  I  )*  M/I  X  0*  I  ) 

V  T  A  =  C  A8S ( VT ) 

IFINl .EQ. 0)G0  TO  9 
VT=  VT  *  2 

VT  A=VT  A*2 

9  WRITEI3.1  12)1  .VII  D.VT.  VTA 

00  11  1=1 .M 

IFINl.  EQ.O)GU  TO  10 
II  =  I  MM 

VT  =  2*CJ*V  1(11) 

V  T  A=C A8  SI VT) 

WRITE  I  3,1  12)1  .VI  (  (I  ) . VT , VTA 
GO  TO  11 
WRITE! 3,1 13)1 


10 
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11  CONTINUE 

WRI TE( 3.1 1 51 
00  13  1=1. MM 
IF  INI .EO. 0) CO  TO  12 
VT=2*CJ*V2! I ) *M/( XO  +  t  ) 
VTA=CABS!  VT) 

WRITE! 3.1 1211 «V2( I  J .VT. VTA 
GO  TO  13 

12  WRITE!  3.1  13)1 

13  CONTINUE 

OO  14  1  =  1  .M 

( 1= (♦MM 
VT=V2!  1  II 
VT  A=C ABS! VT ) 

IF! Nt .EQ. 0 ) GO  TO  14 

VT=VT*2 

VT  A  =  V  T  A  *2 

14  WRITE! 3.1  12II.V2!  II  I.VT.VTA 

STOP 

END 


SAMPLE  OUTPUT 


OUTER  RAOI US/WAVELENGTH  INNER  R AD IU S/ W A VELENGTH  NUMBER  OF 

0.2  000  OE—  0 1  O.OOOOOE  00 

CASE  ANGLE  OF  INCIDENCE  MODES  USED 
IDEGREES) 


1  0.0  1  TO  1 

2  30.0  0  TO  1 


20  POINTS  GAUSSIAN  QUADRATURE  INTEGRATION  IS  USED 


SUBD0MA1 NS 
5 
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EQUIVALENT  MAGNETIC  CURRENT! 

CASE  I 
N*  I 


VERTICAL  POLARIZATION 


1 

0.1  586E-  05 

0.  36  86 E- 02 

0.4758E- 04 

0.11 06E 

00 

0. 

1  l  06  E 

00 

2 

0.2885E-05 

0.6966E-02 

0.4328E-04 

0. 1045E 

00 

0. 

1 0  45  E 

00 

3 

0.4292E-05 

0.  1043E-01 

0.4292E- 04 

0. 1043E 

00 

0. 

1  0  43  E 

00 

A 

0.5626E-05 

0.  1 366 E—  0 1 

0.42  19E-04 

0. 1 025E 

00 

0. 

1025E 

00 

5 

0 • 6880E—  05 

0.  16  71  E—  01 

0.41 28E-04 

0. 1003E 

00 
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0.8829E— 01 

3 

—  0 • 2243E  —  04 

-0.4851  E-Ol 

-0.44  86E-04 

-0.9703E-01 

0.9703E-01 

4 

-0.  2212E-04 

-0. 4655  E-  01 

-0.4424E-  04 

—  0.9310  E—  0 1 

0.9310E-01 

5 

—  0.23  19E  —  04 

—  0.  4764 E—  0 1 

-0.4638E-04 

-  0.9527E-01 

0 • 9527E—  0 1 

6 

—  0.2420E  —  04 

-C.  48 03 E—  0 1 

—  0 . 48  39E—  04 

-0.96 06 E—  0 1 

0.96  06  E—  0 1 

7 

—  0. 2438E—  04 

—  0 . 49 1 7E—  0 1 

-0.48  75E-  04 

-0. 98 33 E- 01 

0. 98  33  E- 01 

3 

-0.2526E-04 

-0.50  1 2  E—  0 1 

—  0  •  50  53  E—  04 

-0.1002E  00 

0.1002E  00 

9 

—  0 • 2684E— 04 

-  0 . 5 1  6  1  E-  0 1 

- 0.5368E-04 

-0.1032E  00 

0.10 32 E  00 

10 

—  0 • 2886E—  04 

—  0  •  53  53  E—  01 

-0.57  72E-04 

-0.1071E  00 

0.1071E  00 

11 

—  0.31  06E—  04 

—  0  •  56  49  E—  0 1 

-0.62  12E-04 

-0.U30E  00 

0.11  30  E  00 

12 

—  0  •  3464E—  04 

—  0 • 6088  E—  01 

—  0.6929E- 04 

—0.121 8E  00 

0.1218E  00 

13 

—  0.4092E—  04 

—  0  •  6887  E—  0 1 

-0.81  84  E—  04 

-0.1377E  00 

0.1377E  00 

14 

-0.4859E-04 

-0.  79  31  E-  01 

- 0.9718E-04 

-0.1586E  00 

0.1586E  00 

15 

-0. 1203E-03 

-0.  1841 E  00 

—  0. 24  05 E—  03 

-0.3682E  00 

0. 3682  E  00 
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Data  for  equivalent  magnetic  current  in  the  above  sample  output 
is  arranged  as  follows.  The  first  and  second  columns  are  the  real  and 
imaginary  parts  of  the  elements  of  vn-L  and  vn//  for  the  two  polarizations. 


For  the  E^-polarization,  the  last  three  columns  are  the  real  part,  the 

m 

e  M  v 

imaginary  part  and  the  magnitude  of  ^'x  an<^  ^3  v^q'  Tor  t^ie  ^  ~ 

polarization,  the  last  three  columns  are  the  real  part,  the  imaginary 

•  vc  n// 

£  j  M  V 

part  and  the  magnitude  of  -.--P*!  and  v”^.  The  significance  of 
these  last  three  columns  is  obvious  from  (117)  to  (119)  and  (132)  to 
(134).  Notice  that  v^1  and  v°#are  always  zero  matrices  and  only  the 


first  two  columns  are  printed  for  them. 
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B-2.1.  Subroutine  FAR 

Subroutine  FAR  (THETA,  PE,  PH,  PEC,  PHC,  TC)  computes  the  following 
two  quantities:  (1)  the  integral 

2  2lT  2 

t(9)  =  rpsln  9  (|F  |2  +  |FQ|2)d0  (B-l) 

in  JQ 

and  (2)  the  ratio 

2 

G'(e,(t>)  =  Hr2  (|Fj2  +  |F0|2)/P.n  (B-2) 

Notice  that  the  transmission  coefficient  and  the  power  gain  can  be  re¬ 
lated  to  the  above  quantities  by 

tt/2 

TC  =  t(6)d0  (B-3) 

and  0 

G(0,<f>)  =  G’(0,<j>)/TC  (B-4) 

Fq,  F,  are  computed  according  to  (147)  to  (150)  and  (151)  to  (154).  The 

u  (p 

integration  in  (B-l)  is  carried  out  analytically  by  integrating  the 
trigonometric  functions. 

This  subroutine  treats  a  number  of  indicences  at  the  same  time. 
t(0)  for  each  incidence  is  computed  and  stored  in  TC.  G'(9,<j>)  for  each 
incidence  is  computed  in  two  different  planes.  If  and  <p^  specify  the 
two  planes  for  the  ith  incidence,  then,  G'(0,<J>^),  G'^,^),  G'(0,4>^+tt) 
and  G' (0 ,(J>2+tt)  are  stored  as  the  ith  element  of  arrays  PE,  PH,  PEC  and  PHC 
respectively.  The  input  variables  are  defined  as 

THETA  =  0,  the  elevation  angle  (in  radians) 

X(i)  =  x^ ,  Y(i)  =  xi_^/,2  as  defined  by  (92) 
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IMS(i)  =  the  order  of  the  lowest  mode  used  for  the 
ith  incidence 

IMF(i)  =  the  order  of  the  highest  mode  used  for  the 
ith  incidence 

PHE(i)  =  angle  specifying  the  first  plane  for  pattern 
of  the  ith  incidence  (in  degrees). 

PHH(i)  =  angle  specifying  the  second  plane  for  pattern  of  the 
ith  incidence  (in  degrees). 

W  stores  the  e2VIU's  starting  with  the  lowest  mode  of 
the  first  incidence. 

Q  stores  the  e^vn//,s  starting  with  the  lowest  mode  of 
the  first  incidence. 

ICASE  =  the  number  of  incidences  considered. 

NMAX  =  the  order  of  the  highest  mode  used  for  all 
incidences  considered. 

M  =  number  of  subdomains  used . 

AK  =  <  as  defined  by  (84). 

RI  =  Rin  inner  radius  (in  wavelengths)  of  the  aperture. 

Minimum  allocations  are  given  by 

COMPLEX  W(NW),  Q(NW),  CFP(NI),  CFT(NI) 

DIMENSION  PE(ICASE),  PH(ICASE),  PEC(ICASE), 

PHC (ICASE) , TC (ICASE) 

COMMON  /BKF1/X  (M-l) ,  Y (M) ,  IMS(ICASE), 

IMF (ICASE) ,  PHE(ICASE) ,  PHH(ICASE) 

/BKF4/BJ (NB) 


mm-*?. 
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CFTn(£) 


j  [J  ...(lac.sin  0)  -  J  ..  (kx  sin  0)] 
n+l  £  n-l  £ 


for  £  = 

"  jn{x£-M+l/2[Jn+l(KX£-M+l/2Sine) + Jn-l(<X £-M+l/2Sln0) ] 
"  *  ^  9  [Jn+2(<y£-W-l/2sln0)  ~ Jn-2(KX £-ll+l/2Sin6)  1 } 


for  £  =  M, . . . ,2M-1 


(B-8) 


DO  loops  5,  6  and  7  compute,  for  the  ith  incidence,  f  (0)  .nd  G ’  ( ©  — 4> ) , 
according  to  the  formulas 

t(0)  =  M(mV&T  IM?(1>  {T  Vn/'12+  cos2e|CFTnTe  vn//|: 

o;  n=IMS(i)  1  i 


+  ~  cos20|CFTn  e2vnI|2  +  | CFP0  e2vni|2} 


(B-9) 


IMF(i)  e 


G’(0,<t>)  =  MfM+7v  ■>  I  N1  cos  n<f>  eivn//+  sin  n<*>  e  vni]T  CFP 

o;  n-IMS(i) 


+  cos20  |  \  '  [sin  n<t>  e. vn//+  cos  ncf>  e-vni]T  CFTn|  } 

n=IMS(l)  1 


(B-10) 


Subroutine  BES,  as  introduced  in  B-1.3  is  called  in  this  subroutine. 
Notice  this  time,  the  output  BJ  of  BES  is  transferred  in  the  common  block 


BLKF4 . 
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LISTING  OF  FAR 


SUBROUT  IN  6  FAR  I  THETA,  PE*  PH.  PEC.  PHC.  TC  ) 

COMMON  /BLKFl /X ( 5 0 > • Y ( 50 ) . t MS ( 2 0 ) . I MF (201 tPHE( 20) . PHH (20) /BLKF2/W. 
6Q/BLKF3/I  CASE  .NMAX*  M.  AK.R  I/BLKF4/BJI  100  ) 

COMPLEX  W (2000) .Q (2 000 ) « CFP (2 50 ).CFT(250) .C J .CK .F TE ,F TH. F PE ,FPH. SP 
6W.SP0. ST* .STQ.FPEC.FPHC.FTEC. FT HC. TERM 
DIMENSION  PEI  20). PH (20 1  .TCI  20 ) . PEC  I  2 0 ) . PHCI 20 ) 

PI=3. 141593 
CJ=(0. . 1. ) 

MM=M—  1 
MD=M»MM 
X0=2*PI*RI/AK 
AK2=AK**2 

AREA=M* (M424X0) /A  K2 
B=AK*S1N(  THETA) 

N1=NMAX+1 
N2— Nl+  1 
N3=N24 1 
BS=B/24 
BC=COS( THETA) 

DO  2  1=1. MM 
T*B*X(  I  ) 

CALL  BESI  N1  ,T  ) 

DO  1  J=1.N  2 

1  BJ I J) =M*B J( J) 

CFTI  I  )  =  2*CJ*BJ(  2) 

CFP  (I  )  =  0. 

K=  l 
CK=l  . 

L  1-2 
L  2=0 

DO  2  J=I,NMAX 
K=K*MD 
L 1 =L 1+1 
L2=L2F 1 
CK=CK*C J 

CFP(K)=CX*IBJ(L1)  +B  J(  L2  )  ) 

2  CFT(K)=CJ*CK* (BJILl )-BJ (L2) ) 

DO  3  1*1.  M 

TX=Y( I ) 

T=B4TX 

CALL  BESIN2.T) 

K=  I  4-MM 

CFPIK  >  =  2*CJ  ♦(  TX*B  J(  2  I  «-BS*(B  J(  D-BJI3)  I  ) 

CFT(K)  =  0. 

X=X4MD 

CFPIK) *-TX4!BJ( 3 ) — B J I l )  )-BS*(3*BJ(2)-BJ (4  )) 

CFT(K)=-C J*(TX*(BJ(3)FBJ(  1)  )-0S • 1 B J I  4 ) F B J 1 2 ) )) 

CK»CJ 
L  1  *3 
L  2*  I 
L  3*4 
L4-2 
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1.5=0 

DO  3  J  =  3.Nl 
CK=CK*C J 
L  l=L  l  4  1 
1.2=1.24  1 
L3=L34  1 
L4=L44 1 
L5=L54 1 
k=*4md 

CFP(KI=CK*CJ*ITX*  (8J(LI  )-Bj(L2  >  >  485* I  2*  0  J  ( L  4  >-B  J <  L 3  )- B J  I L  5 >  >  ) 

3  CFT IK ) =— CK* ( TX* (8 J( LI  )4BJlL2> )— BS*( BJ  <L  3)— B J (L5) ) ) 

K=0 

OO  7  I = 1 • I CASE 
FPE=0. 

FPH=0. 

FTE=0. 

FTH=0. 

FPEC=  0 • 

FPHC=0 . 

FTEC=0. 

FTHC=0 . 

TR  =  0. 

NS=  IMS I 1)41 
NF  =  IMF (1)41 
PH1=PHE(I )*Pt/l 80. 

PH2=PHH ( I >*Pl /I 80 • 

L =MO*  (  NS-  I  ) 

00  6  I  I  =N  S  *  NF 
SP0=0. 

SPW=0. 

STM=0. 

ST  0=0 • 

J=II-t 
SN=C-l  .  )**J 
Sl  =  SIN<  J4PHI ) 

C1  =  C0S<  J* PHI) 

S2=SI N I J4PH2) 

C2=C0S1 J4PH2) 

EN=1. 

F  N=  1  a 

IFIJ.NE.OIGO  TO  4 

EN=.5 

FN=0  • 

4  OO  5  J  J=t  .MO 

L=L4  l 
K*K4l 

TERM=EN4(  Cl  4QIKI4CJ4SI4WIK)  I4CFPILI 
FP  E=FPE4TERM 
FPEC=FP EC 4SN*  TERM 

TERM=EN*(  C2*Q  (K  )4CJ*S2*  W(  K)  )*CFP(L  ) 

FPH=FPH4TERM 
F  PHC=  F  P HC  4  SN*  TE  RM 

T£RM=  EN*I  Cl  *W  IK  >4CJ*S  1  *0<K)  >*CFT  1L  » 

FTE=FTE4TERM 
F TEC=F TEC 4SN* TERM 


TERM=EN*<  C2*W  (K  )*CJ*S  2*000  )*CFT|L) 

fth=fth+t  erm 

F  THC=F  THC  «■  SN*  TERM 
SPO=SPO*a  <K  )*CFP(  L» 

SP*=SPW+W  (K  )*  CFP(  L) 

ST*=  ST  W Fi (K  »*CFT{  L) 

5  STQ=STQ«-(1  (  K  )*CFT<LI 
SPN=EN*CABStSPQ>*  *2  *FN* CABS  CSPW )  **2 
STN=EN*CABS(STW  )♦  *2  *FN*CABS < S TO  1**2 

6  TR=TR  ♦  SPN4-STNA8C*  *2 
FTE=BC*FTE 
FIH®FTH*BC 
FTEC=FTEC*BC 
FTHC=FTHC*BC 

PA  TE=C  ABS  C  F  TE ) *  *2 ♦CAB  S( FPE ) **2 
PATEC=CABS(FTEC  >*  *2 ♦CABS* FPEC ) * *2 
P  A  TH=C ABS  (FTH)**2  *C  ABS ( FPH 1**2 
PATHC  =  CABS<FTHC»*  *2  *-C  AB  S(  FPHC  1**2 
PE ( I ) = PATE *2/ AREA 
PEC ( I )=PATEC*2/AREA 
PH (  I  )=PATH*2/ AREA 
PHC( I )=PATHC*2/AREA 

7  TC(  I  )  =  TR*  SI  N(  THET  A)  /AREA 
RETURN 

END 


B-2.2.  Main  Program  2. 

This  program  computes  the  power  gain  pattern  and  transmission 
coefficients  for  an  aperture.  A  number  of  incidences  are  considered 
at  the  same  time.  Punch  data  cards  are  read  according  to  the  following 
input  sequence: 

101  FORMAT  (13) 

102  FORMAT  (10F6.1) 

103  FORMAT  (2E14.7,  13) 

104  FORMAT  (213) 

105  FORMAT  (5E14.7) 

READ  (1,101)  IT 

READ  (1,101)  ICASE 

READ  (1,102)  (PHE(I) ,  1=1,  ICASE) 

READ  (1,102)  (PHH(I) ,  1=1,  ICASE) 

READ  (1,103)  RO,  RI,  M 

MD  =  2*M-1 
K  =  1 
L  =  MD 

DO  2  1=1,  ICASE 
READ  (1,104)  IMS (I),  IMF(I) 

NSF  =  IMF (I)  -  IMS (I)  +  1 
DO  2  J=l,  NSF 
READ  (1,105)  (W(II),  II=K,L) 

READ  (1,105)  (Q(II) ,  II=K,L) 

K  =  K  +  MD 


2 


L  =  L  +  MD 


The  input  data  is  defined  by 


IT  =  number  of  points  used  in  the  range  of  0  from 
0  to  |  where  pattern  is  computed. 

ICASE  =  number  of  incidences  considered. 

PHE(i)  =  azimuthal  angle  (in  degrees)  specifying  the  first 
plane  where  pattern  is  computed  for  the  ith 
incidence. 

PHH(i)  =  azimuthal  angle  (in  degrees)  specifying  the  second 
plane  where  pattern  is  computed  for  the  ith 
inc idence. 

RI,  RO  are  the  inner  and  outer  radius  (in  wavelengths)  of 
the  aperture. 

M  =  number  of  subdomains  used . 

IMS(i)  =  order  of  the  lowest  mode  used  for  the  ith  incidence. 

IMF(i)  =  order  of  the  highest  mode  used  for  the  ith  incidence. 

W  store  the  e^v111^  starting  from  the  lowest  mode  of 
the  first  incidence 

Q  store  the  e^vn//,s  starting  from  the  lowest  mode 
of  the  first  incidence 
Minimum  allocations  are  given  by 

COMPLEX  W(NW),  Q(NW) 

DIMENSION  PE(ICASE),  PH(ICASE),  PEC(ICASE),  PHC(ICASE), 
TC(ICASE) ,  PTE(NP),  PTH(NP),  TRAN(ICASE) 

COMMON  /BLKF1/X(M-1),  Y(M),  IMS(ICASE),  IMF(ICASE), 


PHE(ICASE),  PHH(ICASE) 
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where 

ICASE 

NW  =  l  (2M-1)  (IMF(i)  -  IMS (i)  +  1) 
i=l 

NP  =  2IT  •  ICASE 

7T 

The  range  0  to  —  of  0  is  divided  into  IT  equal  intervals.  Power  gain 
is  computed  at  the  centers  of  these  intervals  and  the  integration  of 
(B— 3)  is  done  numerically  using  simple  midpoint  rules.  DO  loop  6 
computes  and  store  in  PTE,  PTH  G *  (0 ,4>)  for  2IT  0's  and  both  planes  and 
all  incidences.  The  first  element  of  PTE  is  G'(90°  (1  -  ^)  >  PHE(l))  and 
similarly  for  PTH.  Transmission  coefficients  are  also  computed  in  DO  loop 
6.  TRAN  stores  transmission  coefficients  for  all  incidences.  G(0,$) 
is  computed  according  to  (B-4)  in  DO  loop  8  using  G'(0,<(>)  and  the  trans¬ 


mission  coefficient. 
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LISTING  OF  PROGRAM  2 


COMMON  /aLKFl/X<5  0»  * Y  < 5  0  )  .  I  MS (  2  0  >  .  ( MF (  2 0  I  .PHE!  20 )  ,PHH( 2  0)/8LKF2/W. 
60/BLKF3/I  CASE  ,NMA X, M, AK ,R I 
COMPLEX  W  !  2  00  0)  *0(2000) 

DIMENSION  PE ( 20 )  « PHI  20 )  ,TC( 20 ) , TRAN( 20) . PEC ! 20 ) . PHC ( 2 0) .PTE!  3600 ) . 
6PTHI3600) 

101  FORMAT ( 13) 

102  FORMAT!  10F6.1  ) 

103  FORMAT!2E  14.7.1  3) 

10*  FORMAT  (21  3) 

105  F  OR  MAT  <  5E  1  *  •  7  ) 

106  FORMAT! *1 *. IX, 'FIELD  PATTERN*,//) 

107  FORMAT  1 • 1  », IX ,• CASE* . 13 ./// ,7X,  ‘FIRST  PL ANE •• 3X ,* SECOND  PLANE*./, 
67X.M  •  ,F5.  I  ,‘  DEG.  )  •  ,3X.  •(  *  .  F5 . 1  ,‘DEG.  )•  ./) 

108  FORMAT11X.13.2E1* .♦ ) 

109  FORMAT!//. IX, 'TRANSMISSION  COEFF IC IEN T  = • .E 1 1 • * ) 

READ! 1 . 10 1 ) IT 

READ! 1.101) ICASE 

READ! 1 , 102) !PHE (I >  • 1-1, ICASE) 

READ! 1 • 10  2) (PHH1I  )• 1-1 , (CASE) 

READ! 1 . 10  3) RO  «R ! , M 

PI *3. 1*15 93 

AK*2*P  l*(  RO— R  I  )  /M 

I T 2-=  I  T *2 

MM=M—  1 

MD=M*MM 

D=RO-Rl 

X0*M*R  I  /D 

NMAX=0 

K=  1 

L*MD 

DO  2  1*1. (CASE 

READ!  1.10*)  (MS!  ()  .1  MF  (  I  ) 

NSF=  IMF  !  I  )—  (MS!  0*1 

IF! IMF!  ()  .LE.NMAX >G0  TO  1 

NMAX*t  MF!  I ) 

1  DO  2  J* 1 , NSF 

RE AO! t.l05)!W(((),I I=K,L) 

READ! 1*105) !Q!I() *1 !*K,L) 

K*K*MD 

2  L=L*MD 
T  =  X  0 

DO  3  1*1, MM 
T*T*1 . 

3  X! I )=T 
T*XO— .5 

DO  *  1*1. M 
T*T  *1 . 

*  Y!  I  )*T 

DT*P 1/(2* IT ) 

THETA* ! PI *DT ) /2  « 

DO  5  1*1 . ICASE 
5  TRAN!  0*0. 
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OO  6  1=  It  IT 
THETAs THE TA-OT 

CALL  FAR< THETA.PE .PH.PEC.PHC.TCI 
K=t 

L= I T2* 1 -I 
OO  6  J*l • ICASE 
PTE<KI-PE( Jl 
PTH(K>*PH  {  J  ) 

PTE (L I =PE  C( J» 

PTHIL  )SPHC<  J) 

K*K4I T2 
L=L*IT2 

6  TRANI  JI=TRAN{  J)+DT*TC(J> 

K  =  0 

OO  7  1*1, ICASE 
OO  7  J= I , I T2 

K  =  K-t  I 

PTEIKI=PTE(K)/TRAN< II 

7  PTH(KI*PTH( KI /TRANI  1 1 
MR  I  TE  (3,106) 

K=0 

LL= I- I T  2 

OO  11  I  *1,1  CASE 

WRITE  (3.1  07)1  ,PHE  (  I  ),PHH(  I) 

OO  10  J*1 • I T2 
K=K  +  1 
LL=LL* 1 

10  WRITE(  3.1  08)3. PTE(K).PTHIK) 

11  WRITE ( 3 , 1 09 ) T RANI  I) 

STOP 

ENO 


SAMPLE  OUTPUT 
CASE  1 


FIRST  PLANE  SECONO  PLANE 
I  O.OOEG.)  (  90  •  ODE  G.  I 


1 

0.2195E 

01 

0. 866  OE-O  1 

1  1 

0.218  IE 

01 

0.1 306E 

00 

2 

0.219SE 

01 

0,8  74  IE  — 0  1 

1  2 

0.2178E 

0  1 

0.1 393E 

00 

3 

0,219  4E 

01 

0,8  9  03E**0  1 

13 

0.2175E 

01 

0.1  487E 

00 

♦ 

0.2193E 

01 

0,9  1  45E-0  1 

1  4 

0.217  IE 

01 

0.1 589E 

00 

5 

0.2192E 

01 

0.9468E-0  1 

15 

0.2168E 

01 

0.1697E 

00 

6 

0.2191E 

01 

0,9  870E  — 0 1 

16 

0.2  164E 

01 

0.181  3E 

00 

7 

0,2 18  9E 

01 

0.1  035E  00 

1  7 

0.216  OE 

01 

0.1 936E 

00 

6 

0.2188E 

01 

0,1  09  IE  00 

1  8 

0.2156E 

01 

0.2  066E 

00 

9 

0.2186E 

01 

0.1  1S5E  0  0 

1  9 

0.215IE 

01 
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The  magnetic  current  generated  for  case  2,  E^-polarization, 
in  Section  B-1.4  is  used  to  generate  the  sample  output  above. 
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